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ABSTRACT 


A  numerical  methor  for  calculating  tur¬ 
bulent  axisymmetric  flows:  at  the  stern  and 
in  the  wake  of  bodies  of  revolution  is 
presented.  A  partially  parabolic  marching 
technique  in  a  streamline  coordinate  system 
is  used  together  with  the  k-e  turbulence  model. 

In  the  viscous  sublayer  region,  the  velocity 
is  calculated  using  a  mixing  length  argument 
instead  of  the  wall  function  method.  The 
numerical  procedure  starts  at  a  station  on  the 
body  where  the  boundary  layer  is  thin  and  is 
pursued  several  body  lengths  downstream  into  the 
wake.  The  numerical  solution  is  obtained  by 
marching  downstream  and  iteratively  solving  for 
each  flow  variable.  The  axial  and  normal  cor¬ 
rections  to  the  pressure  are  calculated  by  solving 
a  kinematic  compatibility  equation  for  the  position 
of  the  streamlines.  The  boundary  conditions  for 
the  pressure  are  set  by  calculating  the  potential 
flow  about  an  equivalent  displacement  body.  The 
numerical  marching  scheme  is  repeated,  restarting 
at  the  initial  station,  until  convergence  is 
achieved.  Comparisons  are  made  between  the  numer¬ 
ical  results  and  the  experimental  data  for  four 
different  bodies. 

ADMINISTRATIVE  INFORMATION 

The  work  described  in  this  report  was  funded  by  the  Naval  Sea  Systems  Command 
(NAVSEA  Code  05R24)  Special  Focus  Program  on  Ship  and  Submarine  Drag  and  Wake  Pre¬ 
dictions,  and  was  performed  under  Program  Element  61153N,  Task  Area  SR0230101,  and 
DTNSRDC  Work  Unit  1542-101. 

1 .  INTRODUCTION 

The  present  study  presents  a  numerical  scheme  for  calculating  steady,  axisym¬ 
metric,  turbulent,  incompressible  flow  about  bodies  of  revolution.  To  fully  describe 
the  steady  flow  past  a  body  requires  a  solution  of  the  elliptic  Navier-Stokes 
equations.  Since  this  is  a  difficult  task  for  complicated  turbulent  motion,  the 
turbulence  must  be  modeled.  As  a  minimum  criterion,  the  turbulence  model  must 
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describe  turbulent  motion  effects  on  the  time  averaged  or  mean  flow.  The  present 
numerical  scheme  uses  the  turbulent  kinetic  energy  and  turbulent  kinetic  energy 
dissipation  (k-e)  model  together  with  the  mean  axisymmetric  Navier-Stokes  equations. 

Even  with  a  turbulent  model,  sophisticated  and  time  consuming  numerical  methods 
are  needed  to  solve  the  resulting  elliptical  Navier-Stokes  equations.  However, 
under  certain  conditions  (partially  parabolic  flow  conditions)  certain  diffusive 
terms  are  small.  Neglecting  these  terms,  and  using  a  known  pressure  field,  the 
mean  steady  equations  become  parabolic.  Simpler  marching  techniques  can  then  be 
employed  which  require  less  computer  time  to  calculate  the  solution. 

PARTIALLY  PARABOLIC  FLOW  METHOD 
General  Discussion 

The  partially  parabolic  flow  conditions  are: 

1.  There  is  a  predominant  direction  of  flow; 

2.  The  diffusion  of  momentum,  k  and  e  along  this  direction  can  be 
neglected;  and 

3.  Pressure  is  the  main  carrier  of  downstream  effects  to  the 
upstream  flow. 

A  simple  example  below  illustrates  the  key  points  of  the  partially  parabolic 
flow  assumptions.  For  two-dimensional,  steady,  incompressible,  laminar  flow,  the 
equations  for  the  velocity  components  u  and  v  and  the  normalized  pressure  p  are 
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where  the  first  order  operation  L  is  the  advection  operator.  Assuming  that  the 
partially  parabolic  flow  conditions  are  valid  and  that  the  predominant  direction 
of  the  flow  is  along  the  x-axis,  then  the  momentum  Equations  (1.2)  and  (1.3)  can  be 
approximated  by 


L(u)*-|£+v^  (1.4) 

8y 

L(v)..|£+V0  d.5) 

The  system  of  equations,  given  by  (1.4)  and  (1.5),  is  parabolic  if  the  pressure 
field  is  known.  The  downstream  flow  field  can  then  be  calculated  from  the  u  and  v 
boundary  conditions.  However,  pressure  is  also  unknown;  an  equation  for  p  can  be 
obtained  by  eliminating  the  diffusive  terms  in  Equations  (1.4)  and  (1.5),  and  using 
the  continuity  Equation  (1.1).  The  result  is 


+  2^  -  -  feisi  +  .  (1.6) 

3x2  3y2  L  3x  ^  1 

which  is  a  Poisson  equation.  To  obtain  the  solution  for  the  pressure,  boundary 
conditions  must  be  given  on  the  entire  outer  boundary  of  the  flow  domain,  and  the 
solution  depends  on  upstream  and  downstream  flow  fields.  It  is  directly  through 
the  pressure  solution  that  downstream  effects  are  communicated  to  the  upstream. 

The  partially  parabolic  flow  assumptions  have  been  applied  to  internal  and 
external  flows.  Briley  (1974), 1*  Patankar  et  al.  (1974), 2  and  Roberts  and  Forester 
(1979)^  have  applied  the  assumptions  successfully  to  solve  for  high  Reynolds  number 
flows  in  ducts  and  pipes.  For  external  flows,  the  assumptions  have  been  used  to 
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solve  for  two-dimensional  flows  by  Markatos  and  Wills  (1980),  for  axisymmetric  flow 


*A  complete  listing  of  references  is  given  on  page  ill. 


problems  by  Muraoka  (1980),  and  for  three-dimensional  flow  past  surface  ships  by 
Abdelmeguid  et  al.  (1978)  and  Muraoka  (1978).  For  external  flows,  the  agreement 
with  experimental  data  is  not  as  satisfactory  as  it  is  for  internal  flows.  The 
principle  disagreement  between  calculated  and  experimental  results  for  external 
flows  occurs  in  the  stern/wake  region  of  the  flow.  In  this  region  the  flow  is 
often  near  separation,  and  in  regions  where  the  flow  separates,  the  partially  para¬ 
bolic  flow  equations  cannot  adequately  describe  the  motion. 

Discrepancies  between  experiments  and  numerical  results,  however,  may  also  be 
attributable  to  other  sources.  Among  these  are  the  choice  of  the  predominant 
direction,  the  manner  in  which  the  pressure  is  calculated,  and  the  approximations 
used  in  the  turbulence  model. 

Predominate  Direction  of  Flow 

To  apply  the  partially  parabolic  flow  assumptions,  a  direction  must  be  speci¬ 
fied  along  which  the  diffusion  is  small  compared  with  the  other  terms  in  the 

equations.  Among  the  possibilities  for  this  direction  are  the  axis  of  the  body 
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(Abdelmeguid  et  al.  1978  and  Muraoka  1978,  1980  ),  the  direction  parallel  to  the 

g 

surface  of  the  body  (Pratap  and  Spalding  1975  ),  the  external  potential  flow 

1  4 

direction  (Briley,  1974  ),  and  the  mean  flow  direction  (Markatos  and  Wills,  1980  ). 

With  each  choice  of  the  predominant  direction,  parabolic  equations  are  obtained 

for  momentum,  k  and  e,  but  each  choice  results  in  a  different  approximation  to 

the  turbulent  diffusive  mechanism. 

The  choice  of  the  axis  of  the  body  may  not  be  appropriate  for  bodies  with  a 
steep  turn  at  the  stern.  In  such  a  case,  in  the  stern  region,  the  direction  of  the 
flow  is  sizably  different  from  the  direction  of  the  axis,  and  there  would  be  a  non- 
negligible  contribution  of  diffusion  in  the  axis  direction.  The  selection  of  the 
direction  parallel  to  the  surface  of  the  body,  suitable  for  the  flow  near  the  body, 

becomes  a  poorer  approximation  to  the  flow  direction  as  the  distance  away  from  the 

body  increases.  An  additional  difficulty  with  this  choice  is  that,  at  the  tail  of 
the  body,  the  flow  is  oriented  at  a  nonzero  angle,  while  in  the  wake  the  flow  is 

along  the  axis  of  the  body.  This  can  lead  to  a  singularity  of  the  coordinate  system. 

The  choice  of  the  external  potential  flow  direction  for  external  partially  parabolic 
schemes  is  proper  for  the  external  regions  of  the  flow  but  can  be  a  poor  approxima¬ 
tion  to  the  flow  direction  near  the  body. 


>  W  j  /  V/ 


In  the  present  study,  the  predominant  flow  direction  is  the  mean  flow  direction. 
The  major  advantage  of  this  choice  over  the  others  for  use  in  a  partially  parabolic 
method  is  that  the  direction  of  the  mean  flow  coincides  with  the  same  direction 
that  the  diffusion  is  neglected  in.  As  long  as  the  Reynolds  number  remains  large 
and  the  flow  does  not  separate,  the  diffusion  along  the  mean  flow  direction  is 
small.  Additional  advantages  of  this  selection  are:  there  is  no  singularity  at  the 
tail  of  the  coordinate  system,  and  the  surface  of  the  body  lies  along  a  coordinate 
line.  The  disadvantage  of  this  choice  is  that  additional  equations  are  needed  to 
specify  the  positions  of  the  streamlines. 

Pressure  Field 

The  analysis  leading  to  Equation  (1.6)  shows  that  the  pressure  satisfies  an 

elliptic  differential  equation.  Instead  of  solving  an  equation  of  the  form  (1.6), 

the  partially  parabolic  scheme  employs  the  continuity  equation  to  solve  for  the 

pressure  field.  Returning  to  the  example  of  two-dimensional,  incompressible,  steady, 

*  * 

laminar  flow,  let  u  and  v  represent  the  solution  to  Equations  (1.4)  and  (1.5)  for 

*  *  * 

a  given  pressure  field  p  .  In  general,  u  and  v  will  not  satisfy  the  continuity 

*  *  * 

Equation  (1.1),  and  therefore,  u  ,  v  ,  and  p  must  be  corrected  so  that  the  continu¬ 
ity  equation  is  satisfied.  If  the  corrections  to  the  velocity  and  pressure  are 
denoted  by  primes,  then  the  corrected  u,  v,  and  p  fields  are  given  by: 


u  =  u  +  u 

*  , 
v  =  v  +  v 

* 

P  =  P  +  P 


where  u!  and  v'  must  be  related  by  the  equation 


(1.7) 


5u'  _3v 

3x  3y 


1  /in*  3v^\ 

\3x  3y  / 


(1.8) 


Using  Equations  (1.4)  and  (1.5),  u'  and  v'  can  be  expressed  in  terms  of  the  deriv¬ 
atives  of  p ' : 


-  - '  r) 


(1.9) 


and  inserting  (1.9)  into  (1.8)  leads  to 


(1.10) 


The  functions  f  and  g  in  Equation  (1.9)  must  be  determined  numerically.  The 
solution  to  Equation  (1.10)  together  with  the  forms  of  f  and  g  determine  the  cor¬ 
rected  flow  field. 

Many  partially  parabolic  schemes  have  solved  the  pressure  correction  equation 
by  assuming  that  p'  is  independent  of  the  predominant  direction  of  the  flow.  For 
external  flows,  the  derivatives  of  the  pressure  with  respect  to  the  predominant 
direction  is  set  equal  to  the  value  of  the  derivative  in  this  direction  at  the  edge 
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of  the  boundary  layer  (Patankar  and  Spalding,  1975).  In  this  way,  a  fast  conver¬ 
gent  scheme  can  be  obtained.  However,  treating  of  the  x  and  y  derivatives  of  the 
pressure  as  independent  of  one  another  restricts  the  downstream  influence  of  the 
pressure.  This  can  also  lead  to  poor  results  in  the  thick  boundary  layer  region 
of  the  flow  where  the  longitudinal  and  normal  pressure  variations  can  be  large. 

In  the  present  numerical  scheme,  the  derivatives  of  the  pressure  are  not  assumed 
to  be  independent.  This  slows  down  the  convergence,  but  it  leads  to  a  more  accurate 
calculation  of  the  flow  field.  The  full  details  of  the  present  numerical  scheme 
are  presented  in  the  Numerical  Procedure  Section. 
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For  external  flow  problems,  the  pressure  must  be  specified  on  an  outer 

boundary  away  from  the  body.  Far  away  from  the  surface  of  the  body,  the  pressure 

is  equal  to  the  free-stream  pressure.  This  boundary  condition  has  been  used  by 

4 

Markatos  and  Wills  (1980)  for  two-dimensional  flows.  The  difficulty  with  using 
this  boundary  condition  is  that  special  care  must  be  taken  to  avoid  pressure  oscil¬ 
lations  that  can  develop  in  the  outer  regions  of  the  flow. 

As  an  alternative  to  using  this  boundary  condition,  the  potential  flow  solution 
can  be  used  to  set  the  pressure  on  a  surface  away  from  the  body.  Outside  the 
boundary  layer  region,  viscous  effects  are  negligible  and  the  flow  equations  are 
given  by  the  inviscid  equations.  For  three-dimensional  flows  Abdelmeguid  et  al. 
(1978)6  and  Muraoka  (1978)^  have  applied  the  potential  flow  solution  past  a  surface 
ship  to  set  the  outer  boundary  condition  for  pressure.  In  this  paper,  the  method 
of  using  the  potential  flow  to  set  the  boundary  condition  is  extended  by  using  the 
potential  flow  past  a  displacement  body  (Lighthill,  1958).  The  displacement  body 
concept  is  that  the  actual  body  is  thickened  so  that  outside  the  viscous  flow  region, 
the  streamline  positions  calculated  from  the  potential  flow  about  the  displacement 
body  match  the  streamline  positions  calculated  from  the  viscous  flow  about  the 
actual  body.  Therefore,  the  iteratively  calculated  displacement  body  positions 
and  the  potential  flow  about  this  displacement  body  are  important  calculations  in 
the  numerical  procedure. 

TURBULENCE  MODEL 

The  turbulence  model  used  in  this  study  is  the  k-e  model  (Launder  and  Spalding, 
1974). 11  This  turbulence  model  is  a  two-equation  model,  characterized  by  partial 
differential  equations  for  k  and  for  e.  The  interaction  of  the  turbulence  with 
the  mean  flow  is  specified  through  an  eddy  viscosity  which  is  a  function  of  k  and  e. 

The  k-e  model  has  gained  acceptance  over  mixing  length  models  and  the  k-e  model 
lacks  much  of  the  arbitrariness  of  mixing  length  arguments.  However,  there  is  a 
drawback  to  the  use  of  the  k-e  equations,  that  is:  the  k-e  model  is  not  valid  in 
the  viscous  region  of  the  turbulent  flow.  There  have  been  two  approaches  to 
circumvent  this  difficulty:  (1)  modify  the  k-e  equations  near  the  body  so  that 
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they  are  valid  throughout  the  entire  flow  region  (low  Reynolds  number  model  method, 
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Jones  and  Launder,  1971,  1973  );  and  (2)  instead  of  solving  in  the  viscous  sub¬ 

layer  region,  specify  how  the  velocity  varies  in  this  region  (wall  function  method; 
Launder  and  Spalding  1974'*'''').  The  first  method  is  numerically  time  consuming  and 
the  second  method  does  not  calculate  the  flow  in  a  region  where  the  velocity  has  its 
greatest  variation. 

As  an  alternative  to  these  two  methods,  the  present  study  specifies  the  eddy 
viscosity  near  to  the  body  by  a  mixing  length  model  and  then  calculates  the  velocity 
in  this  region.  In  this  manner  the  problem  of  time  consuming  calculations  is  avoid¬ 
ed  while  the  velocity  is  calculated  in  the  viscous  sublayer  region. 

Another  difficulty  with  solving  the  k-£  model  has  been  found  by  Hanjalic  and 
14 

Launder  (1980  ).  With  partially  parabolic  schemes,  the  k-e  equations  have  been 

solved  by  neglecting  both  the  diffusion  of  k  and  £  in  the  predominant  direction  and 

the  generation  of  k  and  £  due  to  mean  shear  in  the  predominant  direction.  Hanjalic 
14 

and  Launder  (1980  )  retained  the  generation  terms  and  found  considerable  improve¬ 

ment  in  the  consistency  of  the  numerical  results  with  the  experimental  data.  These 
mean  shear  generation  terms  will  be  retained  in  the  k-£  equations  used  in  the 
present  method. 

SUMMARY 

In  summary,  the  present  numerical  scheme  contains  the  following  features  not 
included  in  most  partially  parabolic  schemes  using  the  k-£  turbulence  model. 

1.  A  streamline  coordinate  system  is  used  so  that  the  diffusive  terms  in  the 
flow  equations  are  neglected  in  the  actual  mean  flow  direction; 

2.  The  longitudinal  and  normal  variation  of  the  pressure  field  are  not  treated 
as  independent  quantities; 

3.  The  displacement  body  concept  is  used  to  calculate  the  external  boundary 
condition  for  the  pressure; 

4.  A  mixing  length  argument  is  employed  to  specify  the  eddy  viscosity  in  the 
viscous  sublayer  region;  and 

5.  The  generation  terms  of  k  and  e  due  to  longitudinal  mean  shear  are  retained 
in  the  k-e  equations. 


The  Equation  Section  discusses:  the  geometry  of  the  streamline  coordinate  system, 
the  theoretical  framework  of  the  k-e  turbulence  model,  and  the  time  averaged  equa¬ 
tions.  Also,  the  final  coordinate  system  is  established  together  with  the  equations 
and  their  boundary  conditions.  In  the  Numerical  Procedure  Section,  the  details  of 
the  numerical  scheme  are  given.  Finally,  in  the  Numerical  Results  Section,  results 
for  four  bodies  of  revolution  are  presented  together  with  their  comparisons  with 
experimental  data. 


2.  EQUATIONS 

NATURAL  COORDINATE  SYSTEM 

The  flow  equations  are  expressed  in  the  natural  coordinate  system  so  that  the 
second  partially  parabolic  flow  condition  (see  the  Partially  Parabolic  Flow  Method 
Section)  can  be  applied  in  the  streamline  direction.  A  natural  coordinate  system  is 
one  in  which  one  coordinate  lies  along  the  streamline  and  the  other  two  are  normal 
to  the  streamline  (Liepmann  and  Roshko,  1957). ^  Assuming  axisymmetric  flow,  one  of 
the  normal  directions  is  in  the  azimuthal  direction. 

Differential  distances  in  the  streamline  direction,  the  normal  direction,  and 
the  azimuthal  direction  are  denoted  as  (ds,  dn,  dz) .  If  the  curvilinear  coordinates 
for  this  coordinate  system  are  denoted  by  (£^,  and  the  corresponding  scale 

factors  by  (h^,  r) ,  then  the  differential  distances  are  given  by 

9s  =  hj_  d 

3n  «  h2  d  U  (2.1) 

3z  =  r  d  0  / 


where  8  is  the  azimuthal  angle  and  r  is  the  radial  position, 
any  quantity  are  given  by 


The  derivatives  of 


(2.2) 
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In  the  relation  (2.2),  the  derivative  with  respect  to  0  is  set  equal  to  zero  because 
there  is  no  variation  of  the  flow  in  the  azimuthal  direction  for  axisymmetric  flow. 

In  the  natural  coordinate  system  the  distance  along  the  axis  of  symmetry,  x, 
the  radial  position,  r,  and  the  scale  factors,  h^  and  must  be  determined  in 
terms  of  and  £ As  shown  in  Figure  1  and  using  Equation  (2.1),  the  derivatives 
of  x  and  r  in  the  streamline  direction  are 


=  h^  cos  a 


(2.3a) 


=  h^  sin  a 


(2.3b) 


where  a  is  the  angle  between  the  streamline  and  the  x-axis  direction.  Similarly, 
using  Figure  1  and  Equation  (2.1)  the  derivatives  of  x  and  r  in  the  n  direction  are 
given  by 


-  sin  a 


(2.3c) 


h^  cos  a 


(2.3d) 


An  equation  for  the  scale  factors  can  be  obtained  by  first  differentiating 
Equation  (2.3a)  with  respect  to  and  then  differentiating  Equation  (2.3c)  with 
respect  to  £^.  Performing  these  operations,  the  resulting  left-hand  sides  are  the 
same  and  therefore  the  right-hand  sides  are  equal.  The  result  is: 


(2.4a) 


got  ^^2  Set 

ccs  a  “  hi  sin  a  3^ "  _  iq  sin  a  "  h2  sin  a  si] 


A  second  equation  for  the  scale  factors  can  be  obtained  by  applying  the  same  differ¬ 
entiations  to  Equations  (2.3b)  and  (2.3d),  respectively,  and  the  resulting  equation 


1  .  .  .  Sot 

sin  a  +  h.  cos  a  -57— 

St.2  1  dt*2 


2  ,  .  Sa 

tt=—  cos  a  -  h_  sin  a 


(2.4b) 


A  simpler  equation  can  be  obtained  for  the  scale  factors  by  multiplying 
Equation  (2.4a)  by  -  sin  a  and  multiplying  Equation  (2.4b)  by  cos  a  and  adding  the 
two  resulting  equations,  which  yields: 


3a  _  9h2 

*1  SC, 


(2.5a) 


Similarly,  multiplying  Equation  (2.4a)  by  cos  a  and  multiplying  Equation  (2.4b)  by 
sin  a  and  adding  the  two  resulting  equations  gives  the  result 


3hl  .  3a 
3C2  “  2  3£j_ 


(2.5b) 


Equation  (2.5b)  can  be  put  into  a  more  familiar  form  by  dividing  Equation  (2.5b) 
through  by  h^t^.  Using  Equation  (2.2)  the  result  is 


3a  _  1 
3s  _  R 


(2.6) 


where  R  is  the  radius  of  curvature  of  the  streamlines.  Equations  (2.3)  and  (2.5) 
determine  the  positions  and  the  scale  factors  in  the  natural  coordinate  system. 
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REYNOLDS  EQUATIONS  IN  THE  NATURAL  COORDINATE  SYSTEM  AND 
TURBULENT  EDDY  VISCOSITY 

The  time  averaged,  steady,  axisymmetr ic ,  incompressible  continuity  and  momentum 
equations  will  now  be  given  in  the  natural  coordinate  system  for  the  mean  flow 
fields.  In  this  coordinate  system,  the  variable  s  denotes  the  distance  measured 
along  the  time-averaged  streamline  and  n  denotes  the  variable  measured  normal  to 
this  streamline  direction.  Using  the  results  of  tht.  Natural  Coordinate  System 
Section  the  continuity  and  momentum  equations  are  expressed  as 


h^7  h  (rUh2>  ■  0 


..  3U  1  3p  1  3  ,  — ; — rs  1  9  ,  ,2, 

U  - - -rf- - (r  u'v') - 7T—  (ru'  ) 

9s  p  9s  r  9n  r  9s 


,  ,2  ,2.  1  9r  1  9U  w  9r 

+  (u'  -v'  )  -  77~  +  77  77—  +  -  7— 

r  9s  U  9s  r  9s 


(2.7) 


+  2  u'v’  +  2  v  +  v 
9s  9s2 


9n  \  9n  9s/ 


,  au  A  _2  9r\ 

9s  9n  9sy 

(31  +u  9a\  T  9a  _  1  9rl 

'  9n  9s/  L  9s  r  9nJ 

•5  if 


(2.8) 


1  3£  13 


v'  )  +  w 


2  3s  3n  h.h. 
r  12 


3n  3s  h^r  3n  3s 


In  the  above  equations,  U  is  the  total  time  averaged  velocity  and  its  direction 
coincides  with  the  direction  of  the  mean  streamlines.  The  variable  a  is  the  angle 
between  the  mean  streamline  and  the  direction  of  the  x-axis,  p  is  the  mean  pressure, 
p  is  the  constant  density  of  the  fluid,  V  is  the  molecular  viscosity,  and  r  is  the 
radial  position  from  the  x-axis.  The  velocities  u',  v',  and  w'  are  the  turbulent 
velocities  measured  in  the  s-direction,  in  the  n-direction,  and  in  the  0-direction, 
respectively.  A  bar  above  the  product  of  two  turbulent  velocities  indicates  that 
the  time  average  of  the  product  is  taken.  These  time  averaged  product  terms  are  the 
Reynolds  stress  terms,  which  are  the  stress  on  the  mean  flow  due  to  the  turbulent 
motion.  The  continuity  Equation  (2.7)  states  that  the  mass  contained  between  two 
streamlines,  which  is  r  U  h^  A£^  =  r  U  An,  is  conserved.  Equation  (2.8)  states  that 
the  advection  of  U  in  the  streamline  direction  is  given  by  the  sum  of  the  pressure 
gradient,  the  turbulent  Reynolds  stresses,  and  the  molecular  viscous  forces  in  the 
streamline  direction.  Equation  (2.9)  gives  the  balance  of  the  centrifugal  force 
with  the  normal  pressure  gradient  and  the  normal  turbulent  and  molecular  viscous 
forces . 


To  solve  Equations  (2.8)  and  (2.9),  additional  equations  are  needed  to  describe 
the  Reynolds  stress  terms.  These  terms  are  related  to  the  mean  shear  through  an 


eddy  viscosity  In  cartesian  coordinates,  the  Reynolds  stresses,  -  u|  u^ ,  are 

assumed  to  be  given  by 


u!  u!  +  6  .  ,  k  —  v_ 

i  J  3  ij  T 


(2.10) 


where  the  subscript  on  the  variables  indicates  the  direction  of  the  flow:  (1)  is 
along  the  x-axis,  (2)  is  along  the  y-axis,  and  (3)  is  along  the  z-axis.  Transformed 
to  the  (s,n)  coordinate  system  the  Reynolds  stresses  in  Equation  (2.10)  are  given 
by 


-T—r  3U  LII  3a 

u'v'  =  +U  — 

T  3n  as 


2/3k  -  u’2  =  2VT  § 


~2v„ 


2/ 3k  -  v,Z  =  - -  ~  (r,U) 

r  ds 


2/3K-V'2-  2vtS  || 


-  u ' w'  =  -  v'w'  =  0 


(2.11) 


The  last  two  terms  are  zero  in  Equation  (2.11)  since  they  are  related  to  the  vari¬ 
ation  of  the  mean  flow  in  the  9-direction. 


Combining  Equations  (2.8)  and  (2.9)  with  Equation  (2.11),  the  momentum  equations 


are  written 


U  ™  .  .  i  3*  +  i  3_  L  r  1»1  .  1  L  rv  „  Sal 
3s  p  9s  r  3n  L  e  3n  J  r  3n  [e  9s  J 

h-£]  -(f  £♦*  If)  hff 
-2je“  [f  If ]  2  -  2-e  If  [If  +U  If] 


2Ve  9 [r ,U]*| 
r  9s  J 


(2.12) 


-  i  zt +  7  h  [2ve  h  (rU)] 

1  9  r  9U  JIt  9ct  1  „  U  9r  9r 

-  —  T—  r  v  -T—  +Ur  V  -5—  -  2v  — r  -r—  -r— 

r  3s  e  3n  e  3sJ  e  ^2  3n  3s 

♦ n  [is  ■ ♦  i  h  H 

-*».  [ff«H](f  If) 


(2.13) 


The  variable  v  is  the  total  effective  viscosity,  given  by  the  sum  of  the  molecular 
viscosity  and  the  turbulent  viscosity: 


V  -  V  +  vm 
e  T 


(2.14) 


Once  v^,  is  specified,  a  complete  system  of  equations  is  obtained  for  the  variables 
x,  r,  h^,  h^,  U,  a,  and  p.  The  turbulent  viscosity  will  be  given  by  the  k-e  turbu¬ 
lent  model. 


THE  k-e  TURBULENT  MODEL 

The  k-£  turbulent  model  assumes  that  the  interaction  of  the  turbulence  on  the 
mean  flow  can  be  described  by  an  eddy  viscosity  that  is  a  function  of  the  turbulent 
kinetic  energy  k  and  the  rate  of  turbulent  kinetic  energy  dissipation  e.  In  car¬ 
tesian  coordinates,  k  and  £  are  defined  as 


k  = 


(2.15) 


and 


(2.16) 


where  the  summation  over  the  indices  i  and  j  from  1  to  3  is  implied  in  Equation 

(2.16).  Since  V  is  a  function  of  k  and  £,  it  can  be  shown  by  dimensional  analysis 
A  2 

that  it  is  a  function  of  the  single  variable  k  /£.  It  is  assumed  further  that  the 
function  is  linear  and  vT  is  given  by 

C  k2 

VT  =  (2.17) 

Equations  must  be  supplied  to  describe  the  behavior  of  k  and  e.  The  k  equation 

is  derived  by  considering  the  exact  turbulent  kinetic  energy  budget  which  is  given 
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by  Tennekes  and  Lumley  (1972)  in  cartesian  coordinates  as  ■ 


£ 


(2.18) 


where  summation  is  implied  over  the  indices.  The  left-hand  side  of  (2.18)  is  the 
advection  of  k  by  the  mean  flow.  The  first  three  terms  on  the  right-hand  side  are 
transport  terms:  the  transport  of  k  by  the  turbulent  pressure  gradient,  by  the 
turbulent  velocity  fluctuations,  and  by  the  turbulent  deformation  field.  The  fourth 
term  on  the  right-hand  side  is  the  production  of  k  and  the  fifth  term  is  the  rate 
of  dissipation  of  k.  The  transport  terms  are  parameterized  by  assuming  that  the 
three  bracket  terms  in  Equation  (2.18)  are  proportional  to  the  gradient  of  k: 


—  u'  p'  +  1/2  u'  u!  u!  -  v 
P  j  *  i  j 


ui(V8’‘J 


V, 


T  9k 


ak  3Xj 


(2.19) 


where  is  a  constant.  The  Reynolds  stress  terms  -uj^  uj  in  Equation  (2.18)  are 
given  by  Equation  (2.10).  Combining  these  relations,  the  k  equation  is  given  by 


9k 


j  9x.  9x, 


’^1  H<_1  +  v  +  —L. 

,°k  3*jJ  t[9xj  3xiJ  3xj 


(2.20) 


The  e  equation  is  derived  by  forcing  it  to  take  the  same  form  as  the  k  equation 
(see  Launder  and  Spalding,  1974)^.  It  is  given  by 


9e_  J_  K  3^1  e 

Uj  9xj  ■  9Xj  |^a£  9x^  J  L1  T  k  |^9Xj  9xtJ  9Xj 


C2  e  /k  (2.21) 


where  C^,  C2>  and  a£  are  additional  constants  that  must  be  specified.  The  constants 

C  ,  C. ,  C0,  a,,  and  a  used  to  give  the  best  results  for  boundary  layer  and  wake 
U  1  2  k  e  .. 

.  ,  «  ,,  <  ii.  _  .  j  aOa\  . 


flows  are  given  by  Hanjalic  and  Launder  (1980)  : 


In  the  (s,n)  coordinate  system,  the  k  and  e  equations  are  expressed  as 


and 


3k  _  JL_ 

3_ 

VTh2r  3k 

.  i 

3_ 
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1  W 

|<ro 

_ i 

hir 
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a,  hl  3n 
k 

+  P  -  e  (2.23) 
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+  Cilp-C2F 


(2.24) 


where  the  production  P  of  turbulent  kinetic  energy  is  given  by 


P  = 


+ 


+ 


+ 


2 


? 


(2.2S) 


In  the  viscous  sublayer  region,  Equations  (2.23)  and  (2.24)  do  not  adequately 

describe  the  variations  of  k  and  e.  One  approach  to  remedy  this  problem  is  to  add 

additional  terms  to  Equations  (2.23)  and  (2.24)  to  describe  the  behavior  of  k  and  e 
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near  the  body  (Jones  and  Launder,  1972  ,  1973  ).  This  approach  has  not  been  used 

often  since  the  extra  terms  greatly  increase  the  computational  time  needed  to  cal¬ 
culate  the  flow  field.  Another  approach  has  been- used  to  solve  the  flow  equations 
starting  at  a  point  above  the  viscous  sublayer  region,  using  the  law  of  the  wall  to 
obtain  the  lower  boundary  conditions  for  U,  k,  and  £  (Launder  and  Spalding,  1974).^ 
The  distance  above  the  wall  is  characterized  by 


(2.26a) 


(y-y0)  u, 


(2.26b) 


where  y  is  a  point  above  the  body,  y  is  the  coordinate  of  the  surface,  and  T  is 

o  *  w 

the  shear  stress  at  yQ.  At  a  sufficiently  large  value  of  y+,  given  by  y+  ,  it  is 

assumed  that  the  production  of  turbulent  kinetic  energy  P  is  equal  to  the  rate  of 

dissipation  of  turbulent  kinetic  energy  e.  With  P  =  £,  the  velocity  field  U,  the 

k 

turbulent  kinetic  energy  k,  and  the  rate  of  dissipation  e  are  given  at  y+  =  y+  as 


"  -  “*  [7  +  5-z4] 

k  ■  “*2/^ 


<(y-yQ) 


(2.28a) 


(2.28b) 


(2.28c) 


where  K  is  the  von  Karman  constant  (<»0.40).  Thd  procedure  is  to  determine  <  at 
*  * 
y+  by  solving  Equation  (2.23)  in  the  region  between  y+  *  0  and  y+  =  y+  ,  assuming 

that  the  diffusion  of  k  is  zero.  The  value  of  e  in  Equation  (2.23)  is  set  equal  to 

* 

the  average  energy  dissipation  rate  in  the  region  y  *  0  to  y  =  y  by  integrating 

*  +  +  +  * 
Equation  (2.28c).  Once  k  is  known  at  y+  ,  u^,  U,  and  e  can  be  determined  at  y+ 

from  Equations  (2.28a,  b,  c) .  These  values  of  U,  k,  and  e  serve  as  the  lower 

boundary  conditions  for  the  flow  equations. 

An  alternate  procedure  to  the  two  methods  discussed  above  is  to  specify  the 

* 

viscosity  in  the  region  between  y+  ■  0  to  y+  *  y+  by  a  mixing  length,  eddy  vis¬ 
cosity.  This  is  the  approach  used  in  the  present  study.  A  mixing  length,  eddy 
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viscosity,  denoted  by  v  (whose  form  will  be  given  later),  is  assumed  in  the 

ml  * 

region  between  the  wall  and  y  .  Using  this  eddy  viscosity,  t  can  be  calculated. 

*  +  w  * 

The  value  of  k  at  y+  is  determined  using  Equation  (2.28b),  while  £  at  y+  is  given 

by  Equation  (2.17)  as 

C  k2  * 

e  =  at  y+  =  y+  (2.29) 

Tm£ 

With  these  boundary  conditions,  the  k  and  e  equations  are  solved  in  the  region 
*  * 
y  >  y  .  Using  Equation  (2.29)  as  the  value  of  £  at  y  ensures  that  the  eddy 

T  T  T 

viscosity  calculated  from  the  mixing  length  argument  matches  the  eddy  viscosity 
calculated  from  the  k-£  equations.  This  procedure  provides  a  continuous  variation 
of  the  eddy  viscosity. 


THE  PARTIALLY  PARABOLIC  ASSUMPTIONS  AND  THE  HANJALIC  AND  LAUNDER 
CORRECTION  TO  THE  k-£  MODEL 

The  partially  parabolic  flow  assumptions  given  in  the  Partially  Parabolic  Flow 
Method  Section  will  be  applied  to  the  momentum  Equations  (2.12)  and  (2.13),  and  the 
k-£  Equations  (2.23)  and  (2.24).  The  predominant  direction  in  which  the  diffusion 
terms  are  neglected  is  taken  to  be  along  the  mean  streamline  direction.  Neglecting 
all  s  derivative  terms  in  the  diffusion  terms  of  Equations  (2.12)  and  (2.13),  the 
momentum  equations  are  approximated  by 


1  3£,1  L 

p  3s  r  9n 


(2.30) 


and 


n2  3a  m  _  1  3£ 

9s  p  3n 


(2.31) 
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Patel  et  al.  (1973)  have  shown  that  the  terms  neglected  in  Equations  (2.30)  and 
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(2.31)  are  small  for  axisymmetric  flows.  Using  integral  methods,  Dyne  (1978)  has 
solved  Equations  (2.30)  and  (2.31)  and  has  obtained  good  agreement  with  experimental 
data. 

Equation  (2.30)  states  that  the  momentum  calculated  along  the  streamline 
direction  is  increased  by  the  pressure  force  and  removed  by  the  normal  Reynolds 
stresses.  Equation  (2.31)  states  that  the  centrifugal  force  is  exactly  balanced  by 
the  normal  variations  of  the  pressure. 

The  partially  parabolic  flow  assumptions  are  also  applied  to  the  k  and  e 
equations.  Neglecting  the  turbulent  diffusive  terms  in  Equations  (2.23)  and  (2.24) 
in  the  s-direction,  the  k  and  e  equations  are  approximated  by 


3k 


3s 


hlr 
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3n 


K  1 


3k 

3n 


+  P  -  £ 


(2.32) 


3e 


3s 


V 
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3e 

3n 
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1  k 


P  -  C 


2  k 


(2.33) 


Using  the  partially  parabolic  flow  assumptions  Abdelmeguid  et  al.  (1978)  and 
Muraoka  (1978)^  have  neglected  all  production  terms  containing  derivatives  in  the 


predominant  direction.  In  the  (s,n)  coordinate  system,  with  this  approximation 


P  =  v„ 


m 


(2.34) 
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Hanjalic  and  Launder  (1980)  ,  however,  have  recommended  the  retention  of  the 

(3U/3s)  terms  to  emphasize  the  role  of  irrotational  deformations  in  promoting 
energy  transfer.  Keeping  the  3U/3s  terms  in  P  and  using  Relations  (2.11),  P  is 
approximated  by 


The  difference  in  the  Reynolds  stresses  appearing  in  Equation  (2.35)  is  e: 
pressed  in  terms  of  the  turbulent  kinetic  energy  by  (see  Hanjalic  and  Launder, 
1980)14: 


,  f2  «2\  1  , 

(u'  -v'  )  =  -  k 


(2.36) 


Equations  (2.35)  and  (2.36)  are  used  for  the  production  term  in  Equation  (2.32). 

For  the  e  Equation  (2.33),  Hanjalic  and  Launder  (1980)  4  recommend  the  form  for  the 


production  of  e,  P  ,  as 


P£-ClK(£)2-3  <”'2-V2> 


(2.37) 


where 


C3  =  4.44 


(2.38) 


Combining  Equations  (2.32),  (2.33),  (2.35),  (2.36),  and  (2.37),  the  final  k  and  e 
equations  become 


3k  =  _i_  3  /VT  3k  \  / 3U\  2  1.  30 

3s  h  r  3n  \a,  hl  r  3n  /  VT  \3n)  3  k  9s  "  E 


(2.39) 


“H-V  3^-C2^  (2.40) 


with  the  constants  given  in  Equations  (2.22)  and  (2.38). 
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THE  (x,ip)  COORDINATE  SYSTEM 

The  10  unknown  flow  variables  are  x,  r,  h] ,  h„,  U,  a,  p,  k,  e,  and  v  or  v 

T  mil 

in  the  viscous  sublayer  region.  The  equations  describing  their  behavior  are  (2.3c), 
(2.3d),  (2.5a),  (2.5b),  (2.7),  (2.30),  (2.31),  (2.39),  (2.40),  and  (2.17).  To  inte¬ 
grate  these  equations,  the  curvilinear  coordinates  and  ^  must  be  specified. 

The  coordinate  should  be  chosen  so  that  it  is  the  arc-length  of  some  known  stream¬ 
line.  The  choice  for  can  be  obtained  from  the  continuity  Equation  (2.7).  This 
equation  implies  that  rUh2  is  a  function  of  only  ?2 >  which  is  written 


rUh2 


(2.41) 


The  function  ip  is  the  streamf unction  and  £  should  be  set  equal  to  (p  (Patankar  and 
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Spalding,  1967).  Setting  £2  =  V  a™i  solving  Equation  (2.41)  for  h^  gives  the 
relation 


h 


2 


1_ 

rU 


(2.42) 


For  numerical  calculations,  it  is  advantageous  to  make  one  final  coordinate 
transformation  of  the  equations  in  order  to  replace  by  the  x-axis.  The  trans¬ 
formation  of  the  equations  to  the  (x,ip)  coordinate  system  can  be  obtained  by  con¬ 
sidering  the  variation  of  a  function  f  whose  variables  are  x  and  ip.  The  differential 
of  f  is  given  by 


df 


dx  +  (H) dw 


(2.43) 


Using  Equations  (2.2),  (2.3),  (2.42),  and  (2.43)  the  derivatives  of  f  with  respect 
to  s  and  n  are  given  by: 


3f 


cos  a 


(2.44) 


the  Equations  (for  the  radial  position)  (2.3b)  and  (2.3d),  using  Equations  (2.44) 
and  (2.45),  become 


3r 

-r—  =  tan  a 
3x 


(2.46a) 


and 


3r  =  1 

'd\p  rU  cos  a 


(2.46b) 


Since  ^  is  a  coordinate,  the  conservation  of  mass  is  automatically  satisfied.  In 
place  of  a  conservation  of  mass  equation,  a  kinematic  relation  can  be  obtained  for 
the  (x,y)  coordinate  system  by  differentiating  Equation  (2.46a)  by  y  and  differenti¬ 
ating  Equation  (2.46b)  by  x  and  subtracting  the  two  resulting  equations.  The  result 
is 


3x  (rU  cos  a)  3y  (tan  a)  (2.46c) 

which  expresses  that  the  radial  position  of  the  streamline  calculated  in  the  x- 
direction  from  Equation  (2.46a)  matches  the  radial  position  of  the  streamline  calcu¬ 
lated  in  the  normal  direction  from  Equation  (2.46b).  The  equations  for  U,  a,  k,  and 
e  are  obtained  as  above  by  applying  Equations  (2.44)  and  (2.45)  to  Equations  (2.30), 
(2.31),  (2.39),  and  (2.40).  These  are  given  by 


„  o£  sin  a  3  e  .  3e  sin  a  3 
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(2.46g) 


In  Equations  (2.46f)  and  (2.46g)  the  variations  of  in  the  n-direction  have  been 

neglected  since  they  are  given  by  (3a/3s)  in  Equation  (2.7),  and  this  term  is  taken 

to  be  small  compared  to  the  normal  variations  of  k  and  £.  Outside  the  viscous  sub- 
* 

layer  (y+>y+  )  vT  is  given  by 


vT^cpF 


(2.46h) 


where  y  is  an  intermittency  factor  used  to  reduce  the  turbulent  eddy  viscosity  out¬ 
side  the  boundary  layer  region.  The  factor  is  given  by 


<  0.1 


— - — r  for  t  —  >  0.1 

/r-r  \  6  0  — 


(2.46i) 


1  +  5.5 


where  r  is  the  radial  position  of  the  body  (1  =  0  in  the  wake)  and  6  is  the 

o  o  £ 

boundary  layer  thickness.  In  the  viscous  sublayer  region  (y  <y  )  the  mixing  length 
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eddy  viscosity  is  the  one  used  by  Wang  and  Huang  (1979)  in  the  viscous  sublayer 
region  as  follows: 


where 


l  =  «r 

o 


(2.46k) 


Once  the  boundary  conditions  are  specified,  the  system  of  equations  listed 
under  Equation  (2.46)  are  solved  for  r,  u,  a,  p,  k,  e,  and  vT. 


THE  BOUNDARY  CONDITIONS 

The  solution  to  the  equations  listed  under  (2.46)  is  calculated  in  a  domain  in 
(x,^)  space.  This  domain  begins  at  x^,  which  is  an  x-station  on  the  body  where  the 
boundary  layer  is  thin,  and  ends  at  an  x-station  in  the  wake,  xgl  such  that  xg/ L  = 
4.0,  where  L  is  the  length  of  the  body.  At  xfe  all  flow  variables  must  be  given  so 
that  the  boundary  conditions  at  x  =  xfe  are 


r 

U 
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p  =  p,  }  at  x  =  x. 


(z .47a) 


V-  A  % 


At  xg  the  pressure  is  the  only  variable  that  is  specified  so  that 

p  =  pe  at  x  =  xg  (2.47b) 


The  lower  boundary  is  the  streamline  that  lies  along  the  body's  surface  and  along 
r  =  0  in  the  wake.  This  streamline  is  set  as  =  0.  The  boundary  conditions  on 
the  body  (tpO)  are 


r  =  r 

o 

U  =  0 

a  =  a 
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3n 


Ip 
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3x 
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(2.47c) 


where  rQ  is  the  radial  position  of  the  body's  surface  and  is  the  angle  of  the 

surface  with  respect  to  the  x-axis.  The  pressure  boundary  condition  in  Equation 

(2.47c)  is  obtained  by  setting  U  =  0  in  Equation  (2.46e).  The  boundary  conditions 

for  k  and  e  are  not  given  at  ip  =  0.  These  boundary  conditions  are  given  outside 

the  viscous  sublayer  region  by  Equations  (2.27)  and  (2.29)  on  a  streamline  such 
*  * 

that  y+  =  y+  .  Between  ip  =  0  and  y+  =  y+  the  eddy  viscosity  is  given  by  (2.46j). 
In  the  wake  the  boundary  conditions  on  ip  =  0  are 
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dip 
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iEL  = 
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3k 
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ip  =  0  and 
x/L  >  1 


(2.47d) 


which  express  the  axisymmetric  symmetry  of  the  flow.  For  the  outer  boundary  con¬ 
ditions,  a  streamline  is  selected  such  that  it  lies  entirely  outside  the  turbulent 


boundary  layer  and  turbulent  wake  region  of  the  flow, 
by  ip  and  the  boundary  conditions  on  this  term  are 


k  =  e  =  0 


This  streamline  is  denoted 


(2.47e) 


where  Uq  is  the  far  upstream  constant  velocity  and  pQ  is  the  far  upstream  constant 
pressure. 

The  boundary  conditions  at  x  =  x  will  be  given  from  the  numerical  solution 
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of  the  thin  boundary  layer  equation  obtained  from  the  Wang/Huang  (1979)  boundary 
layer  program.  The  pressure  boundary  conditions  at  x^,  x^,  and  ^  will  be  specified 
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by  a  potential  flow  calculation  (Hess  and  Smith,  1966)  about  the  equivalent  dis¬ 
placement  body.  The  boundary  conditions  specified  by  Equation  (2.47)  are  exactly 
those  needed  for  partially  parabolic  flows,  namely  pressure  given  over  the  entire 
flow  domain  and  U,  k,  and  £  given  on  the  beginning  station  and  on  the  outer  and 
bottom  boundaries. 


3.  NUMERICAL  PROCEDURE 

THE  (x,y)  GRID 

Initially,  the  solution  to  the  boundary  layer  equations  for  the  axisymmetric 

20 

body  is  obtained  from  the  Wang/Huang  (1979)  boundary  layer  program.  This  program 
uses  constant  pressure  profiles,  given  from  a  displacement  body  potential  flow 
calculation,  and  a  mixing  length  eddy  viscosity.  This  viscosity  is  given  by 

v  -  V  for  V  <  vT 
i  mx,  l  o 

(3.1) 

vT  -  f(r  ,6)  0.0168  Ug  6*,  for  vT  >_  vT 
o  i  o 
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where  vT  is  the  inner  viscosity  and  is  the  outer  viscosity.  The  inner  viscosity 

i  o 

is  given  by  Equation  (2.46j).  In  the  outer  viscosity  formula,  U  is  the  potential 

*  ® 

flow  velocity  on  the  displacement  body,  6  is  the  calculated  displacement  body  thick¬ 
ness,  6  is  the  boundary  layer  thickness,  and  f(rQ,<5)  is  a  multiplicative  factor  used 
by  the  program  to  reduce  the  outer  viscosity  in  the  thick  boundary  layer  region, 
given  by 


?o 

The  flow  field  obtained  from  the  Wang/Huang  (1979)  program  is  examined  and  a 
x-station,  x^,  is  selected  where  the  boundary  layer  is  thin  (f(r  ,5)=1).  The  veloc¬ 
ity  profiles  obtained  from  the  calculation  serve  as  the  boundary  conditions  in 
Equation  (2.47a).  The  dissipation  at  x^  is  obtained  by  assuming  P  =  c  and  k  is 
obtained  by  assuming  is  given  by  Equation  (3.1)  at  xfe  and  solving  for  k  from 
Equation  (2.46h). 

The  output  at  x^  is  given  at  J-l  radial  stations.  The  corresponding  stream¬ 
lines  for  this  flow  data  are  calculated  by  integrating  the  definition  of  in  (x,r) 
coordinates: 


3r 


3j£ 

3x 


(3.3) 


setting  the  streamline  on  the  surface  of  the  body  as  =  0.  These  streamlines  are 
denoted  as  (ip ^ ...,  Yj)>  where  ^  =  0.  In  addition  to  the  output  at  x^,  the 
velocity  field  at  x/L  =  1.0  is  examined.  Using  Equation  (3.3),  the  corresponding 
streamlines  are  calculated  at  x/L  =  1.0.  An  outer  streamline  number,  is 

obtained  by  multiplying  1.5  times  the  maximum  of  the  streamline  numbers  obtained 
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at  x^  and  x/L  =  1.0.  The  flow  calculated  on  ip  =  therefore,  should  be  entirely 
above  the  turbulent  region.  The  data  at  x^  is  extended  from  ip  to  by  adding  6 
more  streamline  numbers,  assuming  that  r,  U,  a,  and  p  are  given  by  the  displacement 
body,  potential  flow  solution.  The  streamlines  (y. ,  ....  fv),  where  K  corre- 

sponds  to  the  outer  streamline  number,  serve  as  the  ip  grid  stations.  It  is  empha¬ 
sized  that  these  are  streamline  numbers  and  the  positions  of  these  numbers  must  be 
calculated  from  Equation  (2.46a)  or  (2.46b). 

The  x-stations  are  selected  beginning  at  x^  and  ending  at  xg,  such  that  there 
are  many  x-stations  in  the  stern  and  near-wake  region.  In  this  way,  a  full  resolu¬ 
tion  of  the  thick  boundary  layer  and  wake  region  is  obtained. 

Most  of  the  variables  are  not  evaluated  at  the  grid  points,  but  they  were 
placed  at  staggered  positions  to  give  a  better  representation  for  the  numerical 
calculations.  If  a  typical  grid  point  is  denoted  as  (x1,  ty.),  then  the  flow  vari- 

.  .  .  i+1/2  tT  i+1/2  i  i+1/2  .  i+1/2  i+1/2  .  i+1/2 

ables  are  given  as  ^  ,  Uj+1/2  ,  a.  ,  Pj+1/2  ,  kj+1/2  ,  <-'j+1/2  .  and  vTj+1/2- 

The  i+1/2  superscript  indicates  that  the  variable  is  evaluated  at  (x*+x*+^2)/2  and 
the  j+1/2  subscript  indicates  that  the  position  is  at  (ip.+ip  )/2.  In  addition,  all 


j  rj+l; 


i.  i+1, 


the  variables  are  given  on  ^  a  0  and  As  an  example,  the  U  field  at  (x  +x  )/2 


i+1/2  „  i+1/2  „  i+1/2 


i+1/2  „  i+1/2 


is  given  by  (U/  ....  U^'  ,  UR 

is  illustrated  in  Figure  2. 


.  The  grid  geometry 


THE  FINITE  VOLUME  NUMERICAL  TECHNIQUE 

The  solution  to  the  equations  under  (2.46)  is  obtained  by  marching  the  equations 
downstream  in  x.  Given  the  flow  conditions  at  x^  =  x^;  r,  U,  k,  £,  and  are 
solved  for  at  x^2  aru*  f°r  a  *-s  solved  for  at  x^.  With  these  solutions,  r,  U,  k, 
e,  and  VT  are  obtained  at  x^2  and  a  at  x^  and  so  on  to  the  final  station,  xg.  The 
solution  method  for  the  downstream  values  of  U,  a,  k,  and  £  is  the  finite  volume 
numerical  technique  (Patankar  and  Spalding,  1979). 9  The  values  of  U  at  xi+1/^2  are 
calculated  dividing  the  flow  domain  from  x^  to  x^  ^  into  K-l  areas,  AA  .  ,  and 
integrating  Equation  (2.46d)  in  each  AA^  .  The  control  areas  AA^  are  defined  as 
(see  Figure  2) 


AAj  - 


i-1/2  ^  .  i+1/2 

X  <  X  <  X 


^4  <  V  < 


(3.4) 


VVV\  S  *1 


The  details  of  the  integration  procedure  are  given  in  Appendix  A.  The  result,  at 
each  x-station,  xi+1^2,  is  a  tridiagonal  system  of  equations  for  U (see  Equa¬ 
tion  (A. 6))  which  is  solved  by  the  double  sweep  method  presented  in  Appendix  B.  The 

solutions  for  k  and  £  are  arrived  at  by  the  same  method  except  that  the  lower  control 

* 

areas  begin  at  the  first  streamline  where  y+  >  y+  .  For  smaller  values  of  y+, 

v  =  v  given  by  Equation  (2.46j).  The  value  of  u*  is  calculated  from  the  defini¬ 
ng 

tion  of  T  ,  as  follows: 
w 


2  _  Tw  [  /  TT  9U  .  3u\ 

U*  =  p-  =  [VT  (rU  3n  "  sin  a  3^) 


(3.5) 


V  =  0 


To  solve  for  the  downstream  values  of  a  at  x1,  the  flow  domain  from  x1  1  to  x1 
is  divided  into  K  areas  AA! 1  (see  Figure  2); 


AAj1  = 


i-1  .  .  i 

X  <  X  <  x 


0  £  V  £  (ip^ip^/2 


AA'1  = 


i-1  i 

X  <  X  <  X 


Ip.-l  +  ip.  Ip.  +  \p. 


i-i  .  i 
x  <  x  <  x 


(»WH|’k) 


£  V  £  VK 


for  j  =  2,  ...  K-l 


(3.6) 


and  Equation  (2.46e)  is  integrated  in  each  AA!  .  The  downstream  values  of  r  are 

3  i+1/2 

obtained  by  integrating  Equation  (2.46b)  along  x  =  x  .  The  details  of  the  k, 
e,  a,  and  r  equation  integrations  are  also  given  in  Appendix  A. 


PRESSURE  FIELD  CALCULATIONS 


The  downstream  solutions  Uj^J^  and  a j  ^  depend  on  the  pressure  field. 
Initially,  the  pressure  profiles  are  given  from  the  Wang/Huang  (1979) 20  boundary 
layer  program.  However,  the  downstream  values  of  U  and  a  do  not  satisfy  the  kine¬ 
matic  streamline  relation  in  Equation  (2.46c).  Therefore,  corrections  are  added  to 
the  U,  a,  and  p  values  in  such  a  manner  that  Equation  (2.46c)  is  satisfied.  These 
corrections  are  added  to  the  downstream  values  of  U  and  cl  and  the  upstream  values  of 
p  as  follows: 


U*  i+1/2 
j+1/2 


,  i+1/2 
j+1/2 


+  U 


,  i+1/2 
j+1/2 


\ 


* 

a. 

j 


i 


+  a 


> 


(3.7) 


*  i-1/2  i-1/2  ,  i-1/2 

pj+l/2  "  pj+l/2  +  pj+l/2 


where  the  primed  fields  are  the  corrections  and  the  starred  fields  are  the  new 
corrected  values. 

It  is  the  upstream  values  of  pressure  that  are  corrected  which  is  consistent 
with  the  partially  parabolic  flow  assumption  that  p  communicates  downstream  inform¬ 
ation  to  the  upstream.  Using  the  finite  volume  numerical  technique  to  solve 
Equation  (2.46c),  with  the  control  areas  given  by  Equation  (3.4)  and  with  the  cor¬ 
rected  starred  fields  of  Equation  (3.7),  the  integration  of  Equation  (2.46)  in  the 

AA.  becomes 
J 


,  i+1/2  i-1/2.  r  *  i.  r  *  i 

(x  -x  )  (tan  [oij  ]-tan  ]) 


,  w  r/  i+1/ 2 ,  i+l/2Wo.  *  i+1/2  r  *  i,  *  i  ,o  on 

(qtj -ipj _^) /  t  (rj  +rj-l  Uj-l/2  cos  ^aj  "^j-l^2^  /  (3.8) 


,  w  r  i-1/ 2  i-1/2.,-,  „  i-1/2  r,  i-1  i-1.,-- 

(Y.-<p.  ,)/  [r  +r._.  )/2]  U  cos  t(a.  +a._ ,  )/2. 


The  starred  variables  in  Equation  (3.8)  are  expanded  in  terms  of  Equation  (3.7)  and 
to  first-order  quantities  in  Equation  (3.8)  are  given  by: 


tan  [a.  i]  =  tan  [a  i]  +  a'  i/cos2[a, i] 

3  j  J  3 

cos  [  (cij  *+oij  J  )/2]  *  cos  [  (otj  ^+ctj  *^)/2] 

fi/U*  =  h/tj  ,  ,  i+1/2,.  i+1/2.  , 

U/Uj-l/2  J  U/Uj-l/2  J  j-1/ 2  /CUj-l/2  '  J 


Finally,  and  aj  *  are  related  to  pj  and  Pj-l/2^  by  expanding  the 

solution  equations  for  (A. 6)  and  for  a.^  (A. 10)  in  power  serves  for  U! 

,  i  .  ,  i-1/2  J_i/2  2  J~i/2 

a.  ,  and  Pj_1/2  : 


1+1/2 

3-1/2 


U  ,  i-1/2 , .  U 

3-1/2  P j  —1/2  /(Aj- 


SU-l/2<i/2+Cj-l/2) 


,  i  a  ,  i-1/2  pa  ,  i-1/2 
“j  ~  B3  P3+l/2  +  Cj  Pj-l/2 


(3. 1C 


Inserting  Equations  (3.9)  and  (3.10)  into  (3.8),  a  tridiagonal  system  of  equations 

i-1/2 

are  obtained  for  the  correction  to  the  pressure  field  at  x  : 


P  ,  i-1/2  ,  _P  ,  i-1/2  _P  .  i-1/2  P 

lj-l/2  Pj -3/2  j-1/2  p j  —1/2  +  Cj-l/2  Pj+l/2  j-1/2 


(3.11 


The  boundary  conditions  for  Equation  (3.11)  are  that 


Solving  Equation  (3.11)  for  p’  and  using  Equation  (3.10),  the  corrected  values 
for  U,  a,  and  p  are  obtained. 


NUMERICAL  MARCHING  PROCEDURE 

The  equations  under  (2.46)  are  solved  by  marching  downstream  in  x.  Given  the 
upstream  conditions  at  x1"’1  and  x1-1^2,  the  downstream  flow  is  calculated  at  x 1  and 
xi+l/2.  Once  these  flow  conditions  are  obtained  the  flow  is  calculated  at  xi+1  and 


i+3/2 


below: 


.  i+1/2 

and  so  on  until  x  =  x 


The  numerical  marching  procedure  is  outlined 


1.  Given  the  upstream  values  of  the  flow  at  xi_1  and  xi-1^2,  Ui+1^2  is 
obtained  by  solving  Equation  (A. 6)  and  a1  is  obtained  by  solving 
Equation  (A. 10). 

i-1/2 

2.  The  upstream  values  of  p,  p  ,  are  corrected  by  solving  Equation 
(3.11).  Using  the  new  pressure  field,  the  corrected  downstream  values 
of  U  and  a  are  calculated  from  Equation  (3.10). 

3.  The  downstream  values  of  r,  r*+^  2,  are  calculated  from  Equation  (A. 11) 

4.  The  value  of  u^  is  obtained  from  Equation  (3.5),  V  is  calculated 

m£ 

using  Equation  (2.46j),  k  and  e  are  obtained  by  solving  Equation  (A. 16) 
with  the  boundary  condition  Equations  (2.27)  and  (2.29),  and  VT  is 
given  by  Equation  (2.46h). 

5.  Steps  1  through  4  are  repeated  until  the  flow  fields  converge. 

6.  Steps  1  through  5  are  repeated  at  the  next  downstream  stations,  xi+1 
and  x1+3/2,  until  the  final  x-station  is  reached. 

7.  Beginning  at  x^,  steps  1  through  6  are  repeated  until  all  variables 
converge. 

8.  A  new  displacement  body  is  calculated  and  steps  1  through  7  are 
repeated  until  final  convergence  of  the  flow  field. 

The  maximum  number  of  iterations  at  a  station  (Step  5)  is  set  as  5.  The  new 
displacement  body  (Step  8)  is  calculated  by  matching  the  mass  flow  of  the  turbulent 
flow  from  rQ  to  r^  (the  radial  position  of  y  )  to  the  mass  flow  calculated  from  the 
pressure  field's  potential  flow  from  6  to  r  as  follows: 


pUr  cos  a  dr 


(3.13) 


J 


r 

o 


+  5 


p  U  / 1-c  r  cos  a  dr  = 
o  p 


J 


r 

o 


where 


c 

P 


2 


(P-P0) 


(3.14) 


With  a  new  displacement  body,  new  boundary  conditions  for  p  are  established  by 
recalculating  the  potential  flow  about  the  displacement  body.  The  new  pressure 
field  is  given  as: 


i+1/2 

Pj+l/2 


-  i+1/2 
Pj+l/2 


+  <Pi 


i+1/2  -  i+1/2. 
-pk  > 


(3.15) 


where  p  is  the  pressure  field  from  the  previous  iteration  (old  displacement  body). 
In  correcting  the  upstream  pressure  field  (Step  2)  an  under-relaxation  factor 


has  been  used  to  quicken  the  convergence.  If  p 
-  i-1/2 


i-1/2 

j+1/2 


is  the  uncorrected  pressure 


field  and  p 


j+1/2 


is  the  corrected  pressure  field  obtained  by  solving  Equation 

i-1/2 


(3.11),  then  tie  new  pressure  field  p 


j+1/2 


is  given  by 


P 


i-1/2 

j+1/2 


i-1/2  -  i-1/2 
j+1/2  Pj+l/2 


2 


(3.16) 


With  Equation  (3.16),  it  requires  about  20  total  sweeps  to  obtain  a  solution  that 
converges  to  1%  accuracy. 


NUMERICAL  DIFFICULTIES 

Calculations  were  performed  on  a  Burroughs  B7000:168.  The  number  of  iterations 
needed  to  achieve  a  convergent  solution  depended  on  the  complexity  of  the  body 
geometry.  For  Afterbody  5  (see  the  Numerical  Results  Section),  with  37  x-stations 


and  34  streamline  numbers,  23  iterations  were  required  with  13  minutes  of  processor 
time.  For  Model  A,  however,  only  9  iterations  were  needed  with  6  minutes  of  pro¬ 
cessor  time,  using  the  same  number  of  x  and  \p  stations.  The  region  of  the  flow  field 
where  the  slowest  convergence  occurred  was  the  region  near  the  body  where  the  angle 
of  the  surface  with  respect  to  i|>  is  at  its  maximum.  In  this  region,  the  flow 
rapidly  diverges  away  from  the  body.  As  will  be  seen  in  the  Numerical  Results 
Section,  the  flow  about  Afterbody  5  is  near  separation  while  the  flow  about  Model  A 
is  not.  The  maximum  angle  the  program  could  handle  was  40°  in  the  stern  region. 

Above  this  value,  small  regions  of  reverse  flow  were  seen  to  develop,  as  indicated 
by  the  values  of  the  streamline  angle  a  turning  negative.  After  this  occurred,  no 
convergent  solution  was  obtained.  If  separation  occurs,  it  is  no  longer  advanta¬ 
geous  to  integrate  along  the  streamlines  since  in  the  separated  region  the  values 
of  the  streamlines  become  negative  (assuming  that  the  body  lies  along  the  tfJ=0 
streamline).  In  this  case,  downstream  values  of  the  flow  are  needed  to  obtain  an 
accurate  numerical  solution,  and  the  streamwise  turbulent  diffusion  cannot  be 
neglected  in  comparison  with  the  normal  turbulent  diffusion  terms. 

4 .  NUMERICAL  RESULTS 

Using  the  numerical  procedure  developed  in  Section  3,  calculations  were  per¬ 
formed  for  four  bodies  for  which  experimental  darta  were  available.  These  are 

22 

designated  as  Afterbody  1  (Huang  et  al.,  1979),  Afterbody  5  (Huang  et  al. , 

1980),^  Model  A,  and  Model  B  (Lyon,  1932^  and  1934^).  Table  1  lists  the  following 
geometric  and  flow  parameters  for  each  body:  the  length  of  the  body  L,  the  max¬ 
imum  radius  of  the  body  rmax>  the  upstream  flow  velocity  Uq,  and  the  upstream 
Reynolds  number  Re.  The  aft  body  geometries,  together  with  four  streamline 

i k 

positions,  the  displacement  body,  5  and  the  boundary  layer  thickness,  o,  are 

presented  in  Figures  3  through  6.  In  these  figures,  the  streamline  numbers  are 

2 

normalized  with  respect  to  LU  and  they  correspond  to  (outer)  streamline  grid 

o 

numbers  10,  22,  30,  and  34.  The  aft  body  geometries  for  Afterbody  1  (Figure  3)  and 
Afterbody  5  (Figure  4)  are  characterized  by  parallel  middle  bodies  with  inflected 
sterns.  The  aft  body  geometries  for  Model  A  (Figure  5)  and  Model  B  (Figure  6)  have 
constantly  decreasing  radii,  with  the  curvature  of  the  surface  remaining  convex. 


TABLE  1  -  FLOW  AND  BODY  GEOMETRY  PARAMETERS 


L 

(m) 

r 

max 

(m) 

Uoo 

(m/s) 

R 

e 

Afterbody  1 

3.066 

0.1398 

30.48 

6.60  x  lo6 

Afterbody  5 

2.910 

0.1397 

45.72 

9.30  x  106 

Model  A 

1.778 

0.1778 

17.88 

2.09  x  lo6 

Model  B 

1.778 

0.1778 

17.88 

2.05  x  io6 

The  numerical  results  display  a  thickening  of  the  turbulent  region  in  the  stern/wake 

regions  of  the  flows.  The  displacement  bodies  diverge  significantly  from  the 

physical  bodies  near  the  stern  and  continue  into  the  wake  with  slowly  decreasing 

radii.  In  Figures  3  and  4,  the  computed  displacement  body  and  boundary  layer 

* 

thickness  are  compared  to  the  values  of  6  and  6  obtained  from  the  data  of  Huang 
et  al.,  (1979) ,22  (1980)2"*.  For  both  bodies,  the  computed  5  and  6  lie  slightly 
below  the  experimental  results  in  the  stern/ wake  region,  but  overall  the  agreement 
with  the  experiments  is  good. 

Figures  7  through  28  present  a  detailed  comparison  of  the  computed  flow  field 

* 

to  experimental  results.  In  all  of  these  figures  y+  was  set  equal  to  50.  The 
distributions  of  the  frictional  velocity  u*  and  the  wall-pressure  coefficient  cp  are 
shown  in  Figures  7  through  10  for  the  four  bodies.  The  computed  pressure  distri¬ 
bution  for  Afterbody  1  (Figure  7)  has  a  large  trough  at  the  inflected  stern.  At 
this  region  of  the  body,  the  surface  and  the  streamlines  near  the  surface  have  a 
marked  change  in  curvature.  As  the  streamlines  change  curvature  from  convex  to 
concave,  the  pressure  gradient  changes  from  adverse  to  favorable.  Following  the 
concave  part  of  the  stern,  the  streamline  curvature  becomes  convex  again,  with  a 
corresponding  rise  in  the  pressure  on  the  wall.  The  computed  wall  shear  stress, 
given  by  pu^  ,  drops  rapidly  in  the  adverse  pressure  gradient  region  of  the  flow. 
Accompanying  the  sharp  drop  in  the  wall  pressure,  the  wall  shear  stress  rises 
steeply.  With  the  final  change  of  curvature  of  the  streamlines  at  the  tail,  u*  drop 
dramatically  to  zero.  The  computed  pressure  distribution  agrees  well  with  the 


experimental  data,  having  a  maximum  percentage  difference  of  1%  of  the  total  head 
2 

pu  /2.  However,  the  computed  u.  distribution  reaches  a  smaller  value  than  the 

O  " 

experimental  data,  being  25%  under  the  experimental  result  at  x/L  =  0.97.  (For  the 

pressure  field,  the  percent  difference  between  the  computed  pressure  field,  p 

rcomp 

and  the  experimental  pressure  field,  p  ,  will  be  given  by  100  x  !p  -  d  1/ 

2  exp’  6  }  Kexp  Hcomp 1 

PUQ  •  For  the  velocity  field  the  percentage  difference  between  the  computed 

velocity,  Uco  ,  and  the  experimental  velocity,  Uext>  will  be  given  by  100  x  |Ue 

For  other  variables,  the  percentage  difference  will  be  the  local 

L>  UUii/  U 

difference  between  the  computed  (4>  )  and  the  experimental  (<f>  )  given  by 

comp  exp 

100  x  |<j,  -  «J>  |/|<}>  |  .) 

exp  comp1  1  exp1 

The  computed  ufc  and  c^  distributions  for  Afterbody  5  (Figure  8)  display  the 
same  type  of  behavior  as  exhibited  for  Afterbody  1.  For  Afterbody  5,  the  agreement 
with  the  experimental  data  is  good  for  both  the  wall  frictional  velocity  and  the 
wall  pressure  coefficient.  As  is  evident  by  the  steep  drop  in  u*  in  the  adverse 
pressure  gradient  region  of  the  flow,  the  flow  about  this  body  is  very  near  to 
separation  at  x/L  =  0.93. 

The  computed  u4  and  wall  c  for  flow  past  Model  A  are  given  in  Figure  9.  As 

p 

opposed  to  Afterbodies  1  and  5,  the  stern  of  Model  A  is  not  inflected.  Therefore, 
over  the  aft  region  of  the  body,  the  pressure  gradient  remains  adverse  up  to  the 
tail  of  the  body  and  u^  steadily  decreases  to  zero.  The  computed  wall-pressure 
distribution  lies  slightly  above  the  experimental  results,  with  a  maximum  difference 
between  the  computed  and  experimental  pressure  distributions  of  less  than  3%  at 
x/L  =  0.95. 

The  shape  of  Model  B  is  characterized  by  a  sharply  sloping  stern  region.  The 
values  of  the  surface  angles  near  90°  at  the  stern  caused  some  numerical  difficulty 
for  the  computer  code.  To  obtain  the  computed  flow  field,  only  a  few  x-stations 
were  placed  near  the  tail  of  the  body.  Despite  the  limited  number  of  points,  good 
agreement  of  the  computed  wall  pressure  with  the  experimental  data  is  obtained, 
with  a  maximum  difference  of  less  than  4%  occurring  at  x/L  =  0.95  (Figure  10).  The 
sharp  decrease  in  uA  indicates  that  the  flow  is  nearing  separation  as  the  body 
sharply  turns  downward  at  the  tail. 
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Figures  11  through  16  present  detailed  comparisons  of  the  computed  velocity 
fields  to  the  experimental  results  for  the  four  bodies.  In  Figure  11,  the  computed 
velocity  profiles  show  remarkable  agreement  with  the  experimental  data  for  Afterbody 
1  up  to  the  tail  of  the  body.  The  computed  velocities  near  the  tail  and  in  the 
wake  are  also  in  good  agreement  with  the  experimental  results  (Figure  12).  The 
largest  discrepancy  occurs  immediately  behind  the  body  (x/L=l . 0076) ,  where  the 
computed  axial  velocity  is  8%  less  than  the  experimental  value  at  r  =  0.  The  cause 
of  this  large  difference  is  not  known.  Farther  into  the  wake,  the  agreement  with 
the  experiment  is  very  good,  as  is  evident  by  the  profiles  at  x/L  =  1.1820. 

Computed  and  experimental  velocity  profiles  for  Afterbody  5  are  presented  in 
Figures  13  and  14.  Since  no  measurements  were  made  past  the  body,  comparisons 
cannot  be  made  for  the  wake  velocity  profiles.  For  x-stations  on  the  body,  the 
overall  agreement  with  the  experimental  data  is  good.  At  x/L  =  0.8727  (Figure  13), 
the  computed  radial  velocity  is  5%  under  the  experimental  results,  and  the  computed 
axial  velocity  is  3%  under  the  experimental  data.  Near  the  tail,  at  x/L  =  0.9512 
and  x/L  =  0.9874,  (Figure  14),  the  computed  axial  velocities  are  slightly  greater 
than  the  experimental  profiles  (2%). 

Figures  15  and  16  present  the  total  velocity  profile  results  at  several  x- 
stations  for  the  flows  past  Models  A  and  B.  The  consistency  of  the  computed  velocity 
fields  with  the  experimental  velocities  is  very  good  with  a  maximum  difference  of 
less  than  2%. 

Extensive  measurements  of  the  pressure  field  were  made  for  the  flows  past 
Afterbody  1  (Huang  et  al. ,  1979)^  and  Afterbody  5  (Huang  et  al.,  1980).^  Compar¬ 
isons  of  the  computed  to  the  experimental  pressure  fields  for  the  two  bodies  are 
given  in  Figures  17  through  20.  The  computed  pressure  profiles  and  the  experimental 
pressure  fields  for  Afterbody  1  (Figures  17  and  18)  are  within  1.5%;  for  Afterbody  5 
(Figures  19  and  20)  the  results  differ  by  2%.  At  all  the  x-stations  plotted,  the 
computed  pressure  profiles  agree  very  well  with  the  experimental  results  above  the 
boundary-layer  region.  This  indicates  that  the  displacement  body  concept  has 
correctly  determined  the  influence  of  the  turbulent  flow  on  the  external  inviscid 
flow. 
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The  results  of  the  turbulent  kinetic  energy  calculations  for  Afterbody  1  and 
Afterbody  5,  together  with  experimental  results,  are  given  at  several  x-stations 
in  Figures  21  through  24.  For  x-stations  on  the  bodies,  the  agreement  of  the 
computed  k  with  the  experimental  data  is  good.  These  results  are  encouraging  and 
indicate  the  use  of  the  inner  mixing  length  with  the  k-e  model  gives  a  good 
approximation  to  the  turbulent  field.  In  the  wake,  the  value  of  k  on  the  centerline 
is  its  average  from  r  =  0  to  r  =  0.26.  The  agreement  in  the  wake  for  Afterbody  1 
(Figure  22)  is  poor,  with  the  computed  k  field  25%  under  the  experimental  value  at 
x/L  =  1.1820.  Because  of  this  large  difference,  further  testing  of  the  boundary 
conditions  of  the  k-e  equations  is  needed  in  the  wake. 

The  normal  Reynolds  stress  profiles,  -  u'v',  are  presented  at  several  x-stations 
for  Afterbody  1  and  Afterbody  5  in  Figures  25  through  28.  The  agreement  of  the 
computed  results  to  the  experimental  data  is  fair  in  the  thinner  boundary  layer 
regions  and  becomes  poorer  in  the  thick  boundary  layer/wake  region  of  the  flow. 

For  x-stations  lying  on  the  body,  the  computed  normal  Reynolds  stresses  are  smaller 
than  the  experimental  profiles,  while  in  the  wake  the  computed  values  are  greater 
than  the  experimental  values. 

Figures  29  through  36  present  a  comparison  of  the  partially  parabolic,  k-e 

20 

calculations  to  the  Wang/Huang  mixing  length,  boundary  layer  calculations  for  the 


four  bodies.  Experimental  data  are  also  plotted  in  the  figures.  For  Afterbody  1 
(Figure  29),  the  uA  and  c  distributions  are  in  general  agreement  up  to  x/L  =  0.95. 

20  P 

The  Wang/Huang  calculations  do  not  predict  as  steep  a  drop  in  u*  at  the  concave 


part  of  the  body  as  do  the  k-e  calculations.  In  addition,  there  is  no  trough  at 

20 

the  stern  in  the  Wang/Huang  pressure  distribution.  The  velocity  profiles  computed 


for  Afterbody  1  (Figure  30)  from  the  two  programs  are  consistent  with  the  partially 


parabolic  calculation  agreeing  slightly  better  with  the  experimental  profiles  at 
x/L  =  0.9460.  For  Afterbody  5,  the  partially  parabolic,  k-e  calculations  correctly 


predict  the  steep  drop  in  u^  and  the  pressure  trough  near  the  stern  (Figure  31). 

The  velocity  profile  comparisons  for  Afterbody  5  (Figure  32)  demonstrate  again  the 

20 

agreement  of  the  two  computer  codes  (within  2%),  with  the  Wang/Huang  axial  velocity 
profile  agreeing  slightly  better  with  the  experimental  results  at  x/L  =  0.9874. 
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The  calculations  for  u.  and  c  for  Model  A  (Figure  33)  demonstrate  that  the  calcu- 

P 

lations  agree  up  to  the  tail  of  the  body.  At  the  tail,  the  partially  parabolic 
distribution  drops  rapidly  to  zero  as  opposed  to  the  Wang/Huang-  boundary  layer 
code  which  shows  a  slight  upturn.  Both  pressure  distributions  turn  downward  at  the 
tail,  with  the  partially  parabolic  pressure  distribution  reaching  a  slightly  higher 
maximum.  The  computed  velocity  profiles  for  Model  A  (Figure  34)  are  again  consis¬ 
tent.  The  comparisons  of  uA,  c  ,  and  velocity  for  Model  B  (Figures  35  and  36)  show 
the  same  behavior  as  demonstrated  for  Model  A. 

In  Figures  3  through  36,  y  was  set  equal  to  50.  To  test  the  influence  of 
*  + 

y+  on  the  computed  flow  field,  calculations  were  performed  for  Afterbody  5  with 
* 

y  equal  to  100,  300,  and  500.  The  results  of  these  calculations  are  given  in 

T  * 

Figures  37  through  48.  The  y  =  100  computation  (Figures  37  through  40)  shows 

+  * 

little  change  compared  to  the  y+  =  50  calculation.  Near  the  stern,  there  is  a 
slight  increase  in  u4  (Figure  37).  The  velocity  profiles  (Figure  38)  and  the  pres¬ 
sure  profiles  (Figure  39)  have  not  moved  noticeably,  except  for  the  decrease  in 
pressure  at  x/L  =  0.9874.  However,  the  turbulent  kinetic  energy  profiles  (Figure  40) 

have  marked  dips  near  the  body  at  x-stations  near  the  stern. 

* 

The  uA  distribution,  calculated  for  y+  =  300,  has  further  increased  near  the 
stern  (Figure  41).  In  addition,  the  pressure  at  the  wall  near  the  stern  has  in¬ 
creased.  The  velocity  profiles  (Figure  42)  have  slightly  decreased  from  the  pre- 
* 

vious  case  of  y+  =  100.  This  decrease  in  velocity  is  due  to  the  increase  of 

pressure  near  the  stern,  as  evident  in  Figure  43.  The  turbulent  kinetic  energy  at 

the  stern  (Figure  44)  begins  to  show  erratic  behavior  as  the  matching  point  of  the 

k-£  equations  and  the  mixing  length  is  extended  farther  above  the  body. 

: k 

As  y+  is  further  increased  to  500,  the  computed  uA  distribution  lies  above 

the  experimental  results  (Figure  45).  Both  the  velocity  profiles  (Figure  46)  and 

the  pressure  profiles  (Figure  47)  show  marked  movement  away  from  the  experimental 

* 

values  near  the  stern.  The  erratic  behavior  k  experienced  for  the  case  of  y  =  300 

*  + 
is  even  more  pronounced  for  y  =  500  (Figure  48).  Since  the  object  of  the  study 

+  * 

was  not  to  determine  what  the  optimal  value  of  y  is,  no  definitive  statement  can 

*  + 

be  made  concerning  the  best  value  of  y  to  use.  However,  because  of  the  oscilla- 

+  * 

tory  behavior  of  the  turbulent  fields  as  y  is  increased,  it  is  felt  that  the 
*  +  * 
lower  values  of  y  are  preferable  over  the  higher  values  of  y  . 


Finally,  calculations  were  performed  for  flow  past  Afterbody  5  at  Reynolds 
8  9 

numbers  of  1  x  10  and  1  x  10  .  At  these  two  high  Reynolds  numbers,  the  value  of 
* 

y  required  to  obtain  a  convergent  solution  had  to  be  increased  above  50.  For 
8  * 

Re  -  1  x  10  the  lowest  value  of  y  for  which  a  solution  could  be  obtained  was 

9  +  * 

100,  and  for  Re  =  1  x  10  ,  the  minimum  y+  was  625.  The  calculations  for  these  two 
Reynolds  numbers  were  performed  with  the  same  number  of  x-grid  stations  as  the  case 
of  Re  =  9.3  x  10^,  and  the  number  of  ^-stations  were  increased  to  40.  Comparisons 
of  the  u^  and  c  distributions,  calculated  from  the  partially  parabolic,  k-e  program, 

P  20 

to  the  distributions  calculated  from  the  Wang/Huang  boundary  layer  program  are 

g 

given  in  Figure  49  for  Re  =  1  x  10  .  Near  the  stern,  the  k-e  calculations  predict 
a  36%  smaller  value  of  u^  than  the  mixing  length  theory  (Figure  49).  The  two  pres¬ 
sure  distributions  agree  well  up  to  the  concave  region  of  the  body,  where  the 
partially  parabolic  pressure  distribution  has  a  trough.  The  velocity  profile  com- 

g 

parisons  (Figure  50)  for  Re  =  1  x  10  agree  within  3%.  The  corresponding  compari- 

9 

sons  for  the  case  of  Re  =  1  x  10  are  given  in  Figures  51  and  52.  The  partially 

parabolic,  k-e  code  again  predicts  a  smaller  u^  and  a  pressure  trough  near  the 

stern  (Figure  51).  The  velocity  fields  are  again  in  agreement  (Figure  52),  with  the 

mixing-length  profile  being  5%  smaller  near  the  stern  (Figure  52).  The  values  of 
* 

y  used  for  the  higher  Reynolds  numbers  represent  the  minimum  values  attainable  for 
+  * 
this  particular  code.  The  proper  value  of  y+  as  a  function  of  Reynolds  number 

still  remains  to  be  found.  The  partially  parabolic,  k-e  distributions  at  the  high¬ 
er  Reynolds  numbers  (Figures  49  and  51)  contain  some  roughness.  At  these  higher 

Reynolds  numbers,  more  x-  and  ^-stations  are  needed  to  smooth  these  plots.  However, 

20 

the  number  of  grid  points  were  matched  to  the  Wang/Huang  boundary  layer  program 
to  give  a  true  comparison  of  the  two  computer  codes. 

5.  CONCLUDING  REMARKS 

Overall,  the  agreement  between  the  measured  and  calculated  results  is  encour¬ 
aging.  For  most  of  the  flow  field,  the  velocity,  pressure,  and  k  profiles  are 
correctly  predicted.  The  pressure  and  surface-shear-stress  distributions  also  agree 
well  with  the  experimental  data.  The  major  drawback  to  the  present  scheme  is  that 
the  total  number  of  iterations  needed  for  the  flow  calculations  to  converge  is 
very  large  when  a  large  angle  at  the  stern  is  present. 
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Figure  2  -  The  Staggered  Grid  Geometry 


\p  -  3.05x10 


\p  =  1.30x10“3 


-  5/L 


- 1.17*10^ 

*°**  8x/L 
- 4/  •  2.26x10“< 


Figure  3  -  The  Aft  Body  Geometry,  Four  Streamline  Positions,  the  Displacement  Body  6* 
and  the  Boundary  Layer  Thickness  6  for  Afterbody  1 
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Figure  4  -  The  Aft  Body  Geometry,  Four  Streamline  Positions,  the  Displacement  Body 
Position  6*,  and  the  Boundary  Layer  Thickness  for  Afterbody  5 
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Figure  5  -  The  Aft  Body  Geometry,  Four  Streamline  Positions,  the  Position  of  the 
Displacement  Body  6*,  and  the  Boundary  Layer  Thickness  6  for  Model  A 


Figure  6  -  The  Aft  Body  Geometry,  Four  Streamline  Positions,  the  Position  of  the 
Displacement  Body  <5*,  and  the  Boundary  Layer  Thickness  6  for  Model  B 


The  Plots  for  the  Normalized  Frictional  Velocity  u*  and  the  Pressure 
Distribution  on  the  Body,  c  ,  versus  x/L  for  Afterbody  0 
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The  Velocity  Profiles  at  Three  x-Stations  for  Afterbody 
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gure  16  -  The  Velocity  Profiles  at  Three  x-Stations  for 


k-e  CALCULATIONS 
HUANG,  ET  AL.  (1979) 


urbulent  Kinetic  Energy  Profiles  at  Three  x-Stations  for  After 


Figure  24  -  The  Turbulent  Kinetic  Energy  Profiles  at  Three  x-Stations  for  Afterbody  5 
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Figure  26  -  The  Normal  Reynolds  Stress  Profiles  at  Three  x-Stations  for  Afterbody  1 
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Normal  Reynolds  Stress  Profiles  at  Three  x-Stations  for  Afterbody  5 


CALCI 

IG/HU 


k— e  CALCULATIONS 
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Comparison  of  tf  Total  Velocity  Profiles  at  Two  x-Stations  for  Model  A 


k-6  CALCULATIONS 
WANG/HUANG  CALCULATIONS 
LYON  (1934) 


Comparison  of  the  Total  Velocity  Profiles  at  Two  x-Stations  for  Model  B 
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Figure  38  -  The  Velocity  Profiles  at  Three  x-Stations  for  Afterbody  5  with  y+ 
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Figure  52  -  A  Comparison  of  the  Computed 
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APPENDIX  A 


The  equations  listed  under  (2.46)  are  integrated  using  the  finite  volume 

method.  In  this  approach,  the  values  of  the  downstream  flow  variables  are  obtained 

by  integrating  the  equations  over  small  control  volumes,  as  shown  in  Figure  2,  that 

extend  from  the  bottom  to  the  top  of  the  flow  domain.  For  each  control  volume, 

linear  interpolations  of  the  variables  are  used  to  evaluate  the  integral.  By  way 

of  examples:  if  g  represents  a  field  variable  given  at  the  staggered  grid  points 

i+1/2 


Vj+1/2>  and  denoted  by  gJ+1/2 


(as  U,  p,  k,  £,  v^) ,  and  f  is  a  variable 


given  on  the  grid  points  (x  ,  ^ )  and  denoted  by  fj  (as  a),  then,  for  a  control 

volume  extending  from  xi_^^  to  x1+1/2  antj  fr0m  ^  to  V j  »  three  typical  integrals, 
fg,  f3g/3x,  and  g  3f/3if;  are  evaluated  as: 

i+1/2 

C  V* 

fg  dxdy 


f  P 

P1/2 


(xi+l/2.xi-l/2)  < 


(fP+f,  P) 

)  -J _ l-1-  . 


,  i+1/2,  i-1/2. 

(8j-l/2  +8,j-l/2  ) 


(A.  1) 


and 


1+1/2 

„ x 

I  I 

J  i-1/2 


f  |f  dxdV 


,  \  (fi+f1-l)  ,  i+1/2  i-1/2, 

(Vj-Vj-l)  ~2 -  (gj-l/2  ~8j-l/2  } 


(A. 2) 


i+1/2 


r  p 

J  i-l/2  J  dX 


Vl 


,  i+1/2  i-1/2.  ,  i+1/2.  i-1/2 ...  i  f  i  , 

(x  -x  )  (gj_1/2  +g._1/2  Xfj  -fj.p  (A-3) 
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i-1/2 

For  the  U  momentum  Equation  (2.46),  the  control  volumes  extend  from  x  to 

xi+l/2  an<j  from  y  ^  to  .  Denoting  this  area  by  AA ^  ,  in  Equation  (3.4),  the 
integral  of  each  item  in  this  equation  is  written  as 


I2  +  X3  +  I4  +  X5  +  Z6 


(A. 4) 


In  the  following  integrations,  the  superscript  (n)  on  the  downstream  variables  will 
denote  that  the  variable  is  evaluated  from  the  previous  station  and  that  superscript 
(n+1)  is  the  present  iteration.  For  n  =  0,  we  set  the  zeroth  iterate  equal  to  the 
upstream  value  of  the  variable,  except  for  p,  where  the  value  of  the  initial  guess 
or  the  pressure  result  from  the  previous  total  sweep.  The  factor  (l^^^)*  which 
is  equal  to  1  except  at  i  =  3/2,  where  it  is  zero,  is  multiplied  by  some  of  the 
integral  terms  below  to  indicate  that  these  integrals  are  zero  at  i  =  3/2.  This  is 
done  since  these  integrals  require  the  two  previous  upstream  stations.  The  eval¬ 
uation  of  the  individual  terms  in  Equation  (A. 4)  is 


h  -  jjl 

AA, 


3U 

U  cos  a  tt—  dxdty 
dx 


.  .  r/tl  i+1/2 (n)  i-1/2 Wo1 

<V»j-i>  [<Vi/2  +Vi/2  )/21 


COS 


[aji(n)+  )/2]  (Uj!i/22(n+1)-Uj^i/22)  (A. 5*) 


*2  = 


'I 

AA, 


cos  a  dxdw  = 
dx 


—  (qjj _1>  cos  t  (aj  i^n^j^n^)/2 ; 


,  i+1/2  i-1/2 (n). 

(pj-l/2  "Pj-l/2  } 


(A. 5b) 


T  _  fi  sin  a  3  T  .  .  3u1  .  . 

h  =  JJAa  — —  37  rVe  8in  a  3xJ  dxdV  “ 

f (l-63/2)  sin 

(r j 1+1/2  (n)+r  j  _^1/2  (n>+r j  i-1'2** j ^1/2) /4) 

{[r  i+l/2(n)+r.1+1/2(n))/2]  Veji+J/2(n)  sin  [  (a.  i(n)+aj1{n))/2  ] 

[(Uj-l/22(n+1)‘Uj-l/22)/(xi+1/2“xi"1/2)1 

-t(rji"1/V1-1/2)/2]  v£{'2  sin  [  (0j' 1"1-Kx.f-1)/2] 

[  ^/22-Uj_;?22)/ (x±"1/2-*i"3/2)  ] > 


(A. 5c) 


T  ff  sin  a  3  f  2  .  „  3u  1  ,  , 

I  =  -  - r  V  U  r-  dxc 

4  JJAA  r  3x  L  e  3vJ 


■(1_63/2i){sln[(aji(n)+  ctj^n)/2]/trji+1/2(n)+rj1"1/2+rji_1/2)/A]} 

rUr  i+l/2(n)  1+1/ 2 (n) .  .  ,2  i+l/2(n)„  i+l/2(n) 
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-  (l-63/2i)  [(Uj^22<n)+Uj-l/22)/2]{[(rji+1/2(n>+rji"l/2)/2] 

r/  i+1/2 (n)  i-1/2  i+1/2 (n)^  ,i-l/2.  ,  .  ,  i(n). 

[(vej+l/2  +Ve j  +1  / 2+Ve j -1/2  +Vej-l/2)/4]  sin  (tXj  > 

rfn  i+1/2 (n+1).  i+l/2(n+l)  i-1/2  i-1/2..-, 

l(Uj+l/2  +Uj-l/2  "Uj+l/2  “Uj-l/2  ,/2J 
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(A.5e) 


2[(Uj-l/22<n)+Uj-I/22)/2]{[(rji(n)+r/"1)/2]2 


r,  i+1/2 (n)  i-1/2,  i+1/2 (n)  i-1/2., 

[ (vej+l/2  +vej+l/2+vej-l/2  +vej-l/2)/4] 


i+1/2 (n)  i-1/2  i+1/2 (n)  i-1/2  , 

UUj+l/2  +Uj+l/2  Uj-l/2  +Uj-l/2 


tf„  i+1/2 (n+1)  „  i+1/2 (n+l)./r  “Vj-i 
(Uj+l/2  “Uj-l/2  )/(^j+l  )] 

f,  i(n)  1-1..-.2,,  i+l/2(n)  i-1/2  i+1/2  (n) 

"l(rj-l  +rj-l  )/2]  l(vej-l/2  +vej -1/ 2+vej -3/ 2 


,  i-1/2.  ... 

+vej-3/2)/4 


(A.5f  cont.) 


i+l/2(n)  i-1/2  i+l/2(n)  i-1/2  , 

UUj-l/2  +Uj-l/2  +Uj-3/2  +Uj-3/2  '/4J 


[(Uji-l/2(n+1)^^l/2(n+1))/(^_^_l)} 


(A.5f) 


In  these  integrals,  upstream  differencing  is  used  for  all  derivatives  in  x,  and  1^ 
is  evaluated  so  that  this  term  is  implicit.  In  the  integrals  in  Equation  (A. 5),  the 
y(n+l)  terins  are  conected  and  a  tridiagonal  system  of  equations  results: 


.  U(n+1)  i+l/2(n+l)  _  U(nfl)  lT  i+l/2(n+l) 

j -1/2  Uj -3/2  +  -1/2  j-1/ 2 

✓  .  1  \  /  .  i  \ 


.+1/2 (n+1)  U (n+1 )  n  i+1/2 (n+1) 

■1/2  +  j-1/2  j+1/2 


_  u(n+l)  u (n+1)  .  i+1/2  i-1/2 (n). 

“  Ej-l/2  +  Fj-l/2  (pj-l/2  _Pj-l/2  } 


(A. 6) 


The  coefficients  in  Equation  (A. 6)  depend  on  the  known  upstream  values  and  previous 
itention  values  of  the  downstream  values  of  U,  a,  r,  k,  e,  and  Vg.  Equation  (A. 6) 
is  solved  by  the  double  sweep  method  presented  in  Appendix  B,  given  the  top  and 
bottom  boundary  condition  for  . 

For  the  normal  momentum.  Equation  (2.46e),  the  control  volumes  extend  from 
x*  to  x*,  and  from  i+^j ) / 2  to  (if^+ij^+1) /2 .  The  three  integrals  needed  to  be 
evaluated  in  Equation  (2.46e)  are  written 


J1  *  J2  +  J3 


(A. 7) 


and  if  we  denote  a  control  volume  for  this  equation  by  AA ^ ,  Equation  (3.6),  the 
integral  is 

J1  =11  U2  cos  a  dxdip  *  [  (Vj+i~Vj_i)/2] 

AA , 

J 

[(Uj+l/22+Uj-l/22)/2]2  cos  Kaji“1+aj1(n))/2]  1  (n+1)-aj  1_1)  (A. 8a) 


(A. 8b) 


J 


2 


/  l  i“i\ 
-(x  -x  ) r . 


1-1/ z 


KUj-l/2  +Uj+l/2  )Ui 


,  i-l/2(n)  i-1/2 (n) 

j+1/2  ~Pj-l/2 


) 


and 


J 


3 


sin  a 


l£ 

3x 


dxdy 


[<*J+1-»:1_l)/2]  sin  [(a  i_1-w 1(n,)/2] 
(qj1(n)-q:1i-1(n)) 


(A.8c) 


where 


i(n) 

qi 


(d  i+1/2+n  i+1/2+o  i-l/2(n)  A”1/2 (n)\/A 

(Pj+l/2  +Pj-l/2  +Pj+l/2  +PJ-l/2  )M 


(A.  9) 


The  terms  in  Equation  (A.  8)  are  collected  and  the  equation  for  ot^ 
for  j  *  2,  3,  4....J  as: 


i(n+l) 


is  written 


i(n+l)  a(n+l)  a(n+l)  i-l/2(n)  .  _  a(n+l)  i-1/2 (n) 
'j  +  “j  Pj+l/2  +  Cj  pj -1/2 


The  coefficients  in  Equation  (A. 10)  depend  on  known  quantities,  and,  because  of  the 
nature  of  the  integration,  are  independent  of  the  downstream  values  of  U. 

The  equation  for  the  downstream  radial  position  of  the  streamlines  Equation 
(2.46b)  is  integrated  along  the  x*+^2  line.  But,  first  this  equation  is  multiplied 
by  r  to  obtain  the  form: 


2 

U  cos  a 


Integrating  this  equation  from  Wj  ^  to  ^  along  the  x  yields 


r/+l/2("+l>  -  {(tjW/2<n+l),2  +  2 


i+1/ 2 (n+1) 


cos  [(aji(n)+ctj^n))/2]} 


(A. 11) 


Since  the  r  position  of  ip  =  0  is  known,  the  values  of  all  r ^  are  obtained 

by  summing  Equation  (A. 11). 

The  k  and  e  Equations  (2.46f,  g)  are  solved  simultaneously  as  follows: 


and 


— 

— 

2 

3k  _  sin  a 
3x  r 

9_ 

3x 

rv 

e 

_°k 

,  9k 

sin  a 

9x 

sin  a 

r 

9_ 

9x 

r  v 

e  9k 

ok  9y 

-  U 


3y 


rv 

e  .  9k 
— —  sin  a  •*— 
0,  dx 


■  2 

r  v  ... 
- S.  u 


3U  k2  ^ 

-0.33  C3  e  cos  a|j  +  CyY  f- 


„  9U  .  9U 

rU^-sinaS 


n  2 


-  z 


(A. 12a) 


T, 

U  cos  a  t 
9  x 


sin  a  9 
r  9x 


-vh 


rv 

e  .  9e 
— —  sin  a  -5— 
a  9x 

e 


sin  a  9 
r  9x 


2 

r  v 

- -  U 

a  9y; 


rv 

e  .  9e 
— —  sin  a  -5— 
0  9x 

z 


9U 


♦"If 


r2v  r. 
_ *  9e 

a. 


-0.33  k  cos  a  yk 


„  3U  .  3U 
rU  -5—  -  sin  a  -r— 
9^  ox 


,12 


c2e 


(A. 12b) 


The  control  volumes  for  Equations  (A. 12)  are  the  same  as  those  for  the  U  equation. 
The  integrals  are  written  for  Equations  (A. 12a  and  b)  respectively  as: 


L1  *  L2  +  L3  +  L4  +  L5  +  L6  +  L?  +  Lg  (A. 13a) 

Mx  -  M2  +  Mg  +  M4  +  Mg  +  Mg  +  M?  +  Mg  (A. 13b) 


A  term  in  Equation  (A. 13)  corresponds  to  a  term  in  Equation  (A. 12)  in  the  order  that 
it  is  written.  The  procedure  to  determine  Lg  through  Lg  and  Mg  through  Mg  is  the 
same  as  that  performed  to  determine  lg  through  Ig  in  Equation  (A. 5).  The  only 

difference  is  that  since  the  U,  a,  and  r  equation  integrations  are  performed  before 

.  .  .  .  ,  v  ,  r  „i+l/2(n+l)  i(n+l)  .  i+l/2(n+l) 

the  k  and  e  equation  integrations,  the  values  of  U  ,  a  and  r 

are  known.  These  values  are  used  to  determine  the  coefficients  in  Lg  through  Lg 

and  Mg  through  Mg.  The  v terms,  however,  are  given  from  the  previous  iteration. 

The  last  three  terms  in  Equations  (A. 13a)  and  (A. 13b)  are  determined  by  taking 

2  2 

downstream  values  of  k,  e,  (k  /e)  and  (e  /k) .  To  determine  the  downstream  values 

2  2 

of  the  nonlinear  terms  (k  /e)  and  (e  /k),  a  Taylor  series  in  k  and  e  are  taken  to 
first  order: 

(k2/G)i+1/2  =  (k2/e)i_1/2  +  2(k/e)1-1/2  (ki+1/2-k1_1/2) 

^2^62^1-1/2  (ei+l/2_ei-l/2)  =  2(k/e)i-l/2  ki+l/2 

-  (k2/e2)1-1^2  Gi+1/2  (A. 14a) 


(e2/k) 1+1/2  -  (e2/k) 1-1/2  +  2(e/k) 1-1/2  (ei+1/2-ei-1/2) 

-  (e2/k2)i-1/2  (ki+1/2-k1-1/2)  -  2(e/k)  ci+1/2 

-  r£2/k2)i~1/2  ki+1/2 


(A. 14b) 


With  these  two  approximations,  the  integral  terms  Lg,  L7,  Lg,  Mg,  M7,  and  Mg  are 
given  by: 


L 


6 


0.33  C3 


„  3U  .  . 

e  cos  a  dxdy 


-  0.33  C3 


r  i+1/2 (n+1) 

j-1/2 


cos 


,  i(n+l)  i(n+l)W9l  ,  i+1/ 2 (n+1)  „  i-1/2. 

j  “j-l  CUj-l/2  —  j  —1/2  } 


(A. 15a) 


n  -  JJ  C,y  (k2/ e)  (rU3U/3^-sin  a  8U/3x)2  dxdy 


Av 


(u.  -m  )  (xi+1/2-x1-1/2)  C  V  i+1/2(n+1)r2fk/e)  1-1/2  k  i+1/2<n+1> 

Vj-lMx  *  }  M  Yj-l/2  [2(k/e)j_i/2  kj_1/2 

2.2  i-1/2  i+1/2 (n+1) .  i+l/2(n+l).  i+1/2 (n+1) ... .  i+l/2(n+l) 

-2(k  /e  )j_1/2  ej_1/2  J  U(rj  +rj-i  )/2J  Uj-i/2 

r/TT  i+1/ 2  (n+1)  „  i+1/2  (n+1)  ..  ,  .  ,  .  i(n+lK  i-lWo1 

[(Uj+l/2  -Uj -3/2  /(^j+rVj-l)]  -  8in  [(“j  ^j  )/2] 


/'ll 


i+1/ 2 (n+1) 


(A. 15b) 


L 


8 


edxdy 


-  Opj-^j_1)(xi+1/2-xi-1/2)e 


i+1/2  (..+1) 
j-1/2 


(A. 15c) 


M,  =  -  0.33  C  k  cos  a  dxdy  -  -(V,-V,  ,)  0.33  C 

o  M  dx  j  H 


cos  [ (a 


i(n+l)  i(n+l) 


i+1/2 (n+1)  „  i-1/2, 


j 


+aj-l  )/2UUj-l/2  j-1/2  ) 


i+1/2 (n+1) 
j-1/2 


(A. 15d) 


V  V 


M?  =jj  CL  Y  k(r  U  3U/9y-sin  a  9U/3x)2  dxdy  = 


(W  -V  )(xi+1/2-xi-1/2)  C  C  Y  1+1/2 (n+1)  k  i+l/2(n+l) 
VJ-1MX  X  C1  Cp  Yj-l/2  j-1/2 


{[(r.i+l/2(n+l)  =  rji+1/2(n+l))/2]  „  1 


i+l/2(n+l)  or,TT  i+1/ 2  (n+1) 


3+1/2 


'Uj-3/22<n+1)/%+r'f'j-l)  -  sin  [(ctji(n+1)-Kxji"1)/2]  (Uj^22(n+1) 


-Uj-y22/(xi+1/2-xi"1/2)}2 


(A. 15e) 


M8  “/JAv  C2  e2/k  dxdy  =  -(Wj-^_1)(xi+1/2-xi-i/2)c2 


[2(e/k)j^^2  v^/2(n+i>  -  t*W)£y22  k  *y22(n+L)] 

If  all  the  ki+^-/2(n+l)  an(j  ri+l/2(n+l)  j.pring  are  collected  from  L  through  L„  and 

1  O 

through  Mg,  a  tridiagonal  matrix  system  results,  and  is  given  by: 


(n+1) 

'j-1/2 


i+1/ 2 (n+1)  ‘ 
"j-3/2 


i+1/ 2 (n+1) 


'j-3/2 


»  (n+1) 
Bj  —1/ 2 


i+1/ 2 (n+1) 

‘j-1/2 


i+1/ 2 (n+1) 


j-1/2 


(A. 16) 


,  (n+1) 
Jj  -1/ 2 


i+l/2(n+l) 
J +3/ 2 


i+1/ 2 (n+1) 


'j+3/ 2 


.  (n+J) 
j-1/2 


A  .VvWV  vy*, 


where  Bj  ^1/2^  atM*  Cj-l/2^  are  2x2  matrices>  and  is  a  1  *  2  MtriXi 

Given  Che  Cop  and  boCCom  boundary  condiCions  described  in  Che  CexC,  EquaCion  (A. 16) 
is  solved  using  Che  meChod  presented  in  Appendix  B. 
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APPENDIX  B 


In  the  numerical  integration  procedure,  equations  are  encountered  of  the  form 


B1  -1  +  C1  -2  *  dl 


A.  v.-l  +  B,  v.  +  C.  v . . ,  =  d.,  1  =  2,  3 . N-l 

3  -3  J  “I  3  -3+1  -3  J 


(B.l) 


An  Vjj_„  +  B„  =  d 


1  '  “Nv.t  "N 
— N 


where  A.,  B.,  and  C.  are  n,  x,  n  matrices  and  d.  is  a  column  vector  of  length  n,  all 
3  3  3  “3 

independent  of  the  unknown  column  vectors,  v. 's.  To  solve  Equation  (B.l),  solutions 


-3 


of  the  form  (Carnaham,  Luther,  and  Wilkes,  1969) 


26 


-1 


-i  '  *i  ^i  Ci  — i+1 


(B.2) 


are  sought  where  y^  is  a  column  vector,  is  an  n,  x,  n  matrix  (both  are  to  be 


determined)  and  8^  ^  is  the  matrix  inverse  of  8^-  The  product  8^  is  the  normal 


matrix  multiplication.  Substitution  of  Equation  (B.2)  into  (B.l)  gives  the 


recurrence  relations  for  y^  and  8^ : 


-1 


^i  ~  -  A^  ^i— 1 


(B.3) 


(B.4) 


Rewriting  the  first  equation  of  (B.l)  as 


and  comparing  this  to  the  case  c  -  1  in  Equation  (B.2)  results  in: 


Xl  ■  Bi  1  (B.6) 

Rewriting  the  Nth  equation  of  (B.l)  as 

BN  =  ^  ^N-l  +  \  BN-1  CN-1  Xn 

and  combining  the  vN  terms  yields: 

-N  +  [BN-AN  SN-1  CN-1^  (^N_AN  ^N-l)  =  (B,7) 

which  holds  if  we  set  vN+1  =  0.  Once  the  y^s  and  B/s  are  known,  the  is 
determined  from  Equation  (B.7)  and  vN_^  for  j  X  1  from  Equation  (B.2). 
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NSWC,  White  Oak/Lib 
NSWC,  Dahlgren/Lib 
NUSC/Lib 


Copies 

16  NAVSEA 


1 

SEA 

033 

1 

SEA 

03D 

1 

SEA 

05R 

1 

SEA 

05T 

1 

SEA 

32R 

1 

SEA 

312 

1 

SEA 

55 

1 

SEA 

55N 

1 

SEA 

55W 

1 

SEA 

55W3 

1 

SEA 

56X1 

1 

SEA 

56X4 

1 

SEA 

62P 

3 

SEA 

996 

1  NAVFAC/032C 

1  NADC 


1  NAVENPREDRSCHFAC/ 
Tim  Hogan 


1 

NAVSHIPYD 

PTSMH/Lib 

1 

NAVSHIPYD 

PHILA/Lib 

1 

NAVSHIPYD 

NORVA/Lib 

1 

NAVSHIPYD 

CHASN/Lib 

1 

NAVSHIPYD 

LBEACH/Lib 

2 

NAVSHIPYD 
1  Lib 

1  Code 

MARE 

250 

1 

NAVSHIPYD 

PUGET/Lib 

1 

NAVSHIPYD 

PEARL/Code  202.32 

1  NAVSEC,  NORVA/ 6660 . 03 ,  Blount 

12  DTIC 


1 


1 


AFOSR/NAM 

AFFOL/FYS,  J.  Olsen 


Mr 


MARVD 

1  Div  of  Ship  R&D 
1  Lib 

NASA/HQ/Lib 

NASA/Aiues  Res  Ctr,  Lib 

NASA/Langley  Res  Ctr 
1  Lib 

1  D.  Bushnell 
NBS/Lib 
NSF/Eng  Lib 
LC/Sci  &  Tech 
DOT/Lib  TAD-491.1 
MMA 

1  National  Maritime  Res  Ctr 
1  Lib 

U  of  Bridgeport/E.  Uram 
Brown  Univ/J.T.C.  Liu 


U  of  Cal.,  San  Diego 
1  A.T.  Ellis 
1  Scripps  Inst  Lib 

CIT 

1  Aero  Lib 
1  T.Y.  Wu 
1  A. J.  Acosta 
1  I.  Saber sky 
1  D.  Coles 

City  College,  Wave  Hill/Pierson 

Catholic  U  of  Amer/Civil  & 

Mech  Eng 

Colorado  State  U/Eng  Res  Ctr 


1  U  of  Connecticut/Scottron 

1  Cornell  U/Shen 

1  Florida  Atlantic  U/Tech  Lib 

2  Harvard  U 

1  G.  Carrier 

1  Gordon  McKay  Lib 

1  U  of  Hawaii/Bretschneider 

1  U  of  Illinois/ J.  Robertson 

4  U  of  Iowa 

1  Lib 

1  L .  Landweber 
1  J .  Kennedy 
1  V.C.  Patel 

1  Johns  Hopkins  U/Lib 

1  Kansas  State  U/Nesmith 

1  U  of  Kansas/Civil  Eng  Lib 

1  Lehigh  U/Fritz  Eng  Lab  Lib 

MIT 

1  Lib 

1  J . R .  Kerwin 
1  P .  Leehey 

1  J.N.  Newman 

2  U  of  Minn/St.  Anthony  Falls 

1  Lib 
1  R.  Arndt 

1  U  of  Mich/NAME/Lib 

1  U  of  Notre  Dame/ Eng  Lib 

1  New  York  U/Courant  Inst/Lib 

4  Penn  State 

1  B.R.  Parkin 
1  R.E.  Henderson 
1  ARL  Lib 
1  G.H.  Hoffman 


U  of  Cal/Dept  Naval  Arch,  Berkeley  4 
1  Lib 

1  W.  Webster 


Copies 


Copies 


1 

1 

1 

1 

1 

2 


Princeton  U/Mellor 

U  of  Rhode  Island/F.M.  White 

SIT/Lib 

U  of  Texas/Arl  Lib 
Utah  State  U/Jeppson 


Southwest  Res  Inst 
1  Applied  Mech  Rev 
1  Abramson 


Stanford  U 
1  Eng  Lib 

1  R.  Street,  Dept  Civil  Eng 
1  S.J.  Kline,  Dept  Mech  Eng 


Stanford  Res  Inst/Lib 
U  of  Virginia/Aero  Eng  Dept 
U  of  Washington/Arl  Tech  Lib 


VP  I 

1  Dept  Mech  Eng 
1  J.  Schetz,  Dept  Aero  & 
Ocean  Eng 


Webb  Inst 
1  Lib 
1  Ward 


Woods  Hole/Ocean  Eng 
Worchester  Pi/Tech  Lib 
SNAME/Tech  Lib 
Bell  Aerospace 

Bethlehem  Steel/Sparrows  Point 
Bethlehem  Steel/New  York/Lib 
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1 

1 

1 

1 

1 

1 

1 

4 


1 

1 

1 

1 


1 

1 

1 

1 

1 


Boeing  Company/Seattle 
1  Marine  System 
1  P.  Rubber t 


Bolt,  Beranek  &  Newman/Lib 
Exxon,  NY/Design  Div/Tank  Dept 
Exxon  Math  &  System,  Inc . 
General  Dynamics,  EB/Boatwright 
Flow  Research 
Gibbs  &  Cox/Tech  Info 
Grumman  Aerospace  Corp/Lib 


Hydronautics 
1  Lib 

1  E.  Miller 
1  V.  Johnson 
1  C . C .  Hsu 


Lockheed,  Sunnyvale/Wald 
Lockheed,  California/Lib 
Lockheed,  Georgia/Lib 


McDonnell  Douglas,  Long  Beach 
1  T.  Cebeci 


Newport  News  Shipbuilding/Lib 
Nielsen  Eng  &  Res 
Northrop  Corp/Aircraf t  Div 
Rand  Corp 


Rockwell  International 
1  B.  Ujihara 


SAI  (C.  von  Kerczek) 
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Copies 

1 

1 

1 

1 

1 

1 

1 

1 


Copies 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

30 

3 

1 

1 

1 


Sperry  Rand/Tech  Lib 
SRA  (S.J.  Shamroth) 

Sun  Shipbuilding/Chief  Naval  Arch 

Robert  Taggert 

TRW  Systems  Group/Lib 

TRACOR 

United  Techno logy /East  Hartford, 
Conn 

Westinghouse  Electric 
1  Gulino 


CENTER  DISTRIBUTION 


Code 

Name 

1500 

W.B.  Morgan 

1504 

V.J.  Monacella 

1507 

D.S.  Cieslowski 

1508 

R.S.  Rothblum 

152 

W.C.  Lin 

1521 

P.  Pien 

1521 

W.  Day 

1522 

G.  Do bay 

1522 

M.  Wilson 

154 

J.  McCarthy 

154 

P.  Granville 

1540.1 

B.  Yim 

1540.2 

R.  Cumming 

1542 

T.T.  Huang 

1542 

N.C.  Groves 

1542 

G.S.  Belt 

1542 

Y.T.  Lee 

1542 

M.S.  Chang 

Copies 

Code 

Name 

1 

1544 

S.  Jessup 

1 

1544 

K.F.  Lin 

1 

1560 

M.  Martin 

1 

1561 

G .  Cox 

1 

1563 

W.E.  Smith 

1 

1564 

J.  Feldman 

1 

1572 

E.  Zarnick 

1 

1606 

T.C.  Tai 

1 

1615 

R.J.  Furey 

1 

1802.1 

H.  Lugt 

1 

1802.2 

F.  Frenkiel 

1 

1840 

J.  Schot 

1 

1843 

H.  Haussling 

1 

1843 

K.  Rao 

1 

19 

M.M.  Sevik 

1 

1940 

J . T .  Shen 

1 

1942 

B.E.  Bowers 

1 

1942 

T.M.  Farabee 

1 

1942 

F.E.  Geib 

1 

1942 

T.C.  Mathews 

10 

5211.1 

Reports  Distribution 

1 

522.1 

Unclassified  Lib 

(C) 

1 

522.2 

Unclassified  Lib 

(A) 

1 

1 

1 


1544 

1544 

1544 


T.  Brockett 
R.  Boswell 
E.  Caster 


